如何在scipy中使用离散变量值最小化函数

Question

我正在尝试优化具有多个输入变量(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()

Answer 1

你试图做的（如果我理解你的设置）被称为整数规划，它是 NP-hard；http://en.wikipedia.org/wiki/Integer_programming。我意识到您不是在寻找整数解，但是如果您将所有输入乘以 10 并将目标函数除以 100，则会得到一个等效的问题，其中输入都是整数。关键是，您的输入是离散的。

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

有一些用于整数规划问题的近似算法，您可能会发现有时其中一种算法效果很好，可以让您达到最佳点。在我引用的维基百科文章中列出了一些您可以尝试的事情，但我不知道您是否会乐于尝试解决这个问题。