9

我正在尝试优化具有多个输入变量（24 到 30 之间）的目标函数。这些变量是三个不同统计变量的样本，目标函数值是t检验概率值。误差函数表示期望和实际 t 检验概率之间的误差（差的平方和）。对于所有三个 t 检验，我只能接受误差小于 1e-8 的解决方案。

我正在使用scipy.optimize.fmin，效果很好。目标函数为零的解决方案有很多。

问题是我需要找到一个解决方案，其中变量在 0 到 10.0 之间，并且是整数或不超过一位小数部分。有效值的示例是0 10 3 5.5 6.8。无效值的示例：-3 2.23 30或0.16666667。

我碰巧知道至少有一种解决方案，因为目标值来自实际测量数据。原始数据丢失了，我的任务是找到它们。但我不知道怎么做。使用试验/错误不是一种选择，因为每个变量大约有 100 个可能的值，并且考虑到变量的数量，可能的情况数将是 100**30，这太多了。使用 fmin 很棒，但是它不适用于谨慎的值。

有没有办法解决这个问题？如果我需要运行一个程序数小时才能找到解决方案，这不是问题。但是我需要在几天内找到大约 10 个目标值的解决方案，而且我没有新的想法。

这是一个示例 MWE：

import math
import numpy
import scipy.optimize
import scipy.stats
import sys

def log(s):
    sys.stdout.write(str(s))
    sys.stdout.flush()

# List of target T values: TAB, TCA, TCB
TARGETS = numpy.array([
    [0.05456834,   0.01510358,    0.15223353   ],  # task 1 to solve
    [0.15891875,   0.0083665,     0.00040262   ],  # task 2 to solve
])
MAX_ERR = 1e-10 # Maximum error in T values
NMIN,NMAX = 8,10 # Number of samples for T probes. Inclusive.

def fsq(x, t, n):
    """Returns the differences between the target and the actual values."""
    a,b,c = x[0:n],x[n:2*n],x[2*n:3*n]
    results = numpy.array([
        scipy.stats.ttest_rel(a,b)[1], # ab
        scipy.stats.ttest_rel(c,a)[1], # ca
        scipy.stats.ttest_rel(c,b)[1]  # cb
    ])
    # Sum of squares of diffs
    return (results - t)

def f(x, t, n):
    """This is the target function that needs to be minimized."""
    return (fsq(x,t,n)**2).sum()

def main():
    for tidx,t in enumerate(TARGETS):
        print "============================================="
        print "Target %d/%d"%(tidx+1,len(TARGETS))
        for n in range(NMIN,NMAX+1):
            log(" => n=%s "%n)
            successful = False
            tries = 0
            factor = 0.1
            while not successful:
                x0 = numpy.random.random(3*n) * factor
                x = scipy.optimize.fmin(f,x0, [t,n], xtol=MAX_ERR, ftol=MAX_ERR )
                diffs = fsq(x,t,n)
                successful = (numpy.abs(diffs)<MAX_ERR).all()
                if successful:
                    log(" OK, error=[%s,%s,%s]\n"%(diffs[0],diffs[1],diffs[2]))
                    print " SOLUTION FOUND "
                    print x
                else:
                    tries += 1
                    log(" FAILED, tries=%d\n"%tries)
                    print diffs
                    factor += 0.1
                    if tries>5:
                        print "!!!!!!!!!!!! GIVING UP !!!!!!!!!!!"
                        break
if __name__ == "__main__":
    main()
4

1 回答 1

6

您正在尝试做的事情（如果我理解您的设置）称为整数编程，它是 NP 难的；http://en.wikipedia.org/wiki/Integer_programming。我意识到您不是在寻找整数解决方案，但是如果您将所有输入乘以 10 并将目标函数除以 100，您会得到一个输入都是整数的等效问题。关键是，您的输入是离散的。

您正在使用的目标函数是一个凸二次函数，并且有很好的约束优化算法可以快速求解区间 [0, 10] 中的实值输入。从此您可以尝试舍入或检查附近所有可接受的点，但其中有 2^n 个，其中 n 是输入的数量。即使您这样做，也不能保证最佳解决方案是这些点之一。

整数规划问题有一些近似算法，您可能会发现有时其中一种算法运行良好，足以让您达到最佳点。我引用的 Wikipedia 文章中列出了您可以尝试的方法，但我不知道您是否会乐于尝试解决此问题。

于 2013-09-13T05:28:16.400 回答