Scipy 最小化整个数据集和除零以外的所有值的约束答案

【问题标题】：Scipy minimize with constrain for whole data set and all values except zeroScipy 最小化整个数据集和除零以外的所有值的约束
【发布时间】：2021-03-02 19:04:55
【问题描述】：

我尝试通过 most log-likelihoode 方法为我的数据获取统计方程的参数。我想适应约束的方程如下所示： Equation and Constraints

最佳拟合是在我的似然函数的最小值处实现的l

一个最小的示例如下所示：

import numpy as np

from numpy import log
from scipy.optimize import minimize

y =np.array([0.21863326, 0.19805154, 0.22953017, 0.21906749, 0.22067327,
       0.20638931, 0.17845443, 0.20008429, 0.21702199, 0.16334912,
       0.18480577, 0.17172182, 0.16495525, 0.15907978, 0.17029919,
       0.14020628, 0.16528562, 0.17141436, 0.14978351, 0.1329871 ,
       0.14036109, 0.15894933, 0.16783223, 0.17372222, 0.15986161,
       0.1654368 , 0.16348146, 0.15595923, 0.15192792, 0.12272897,
       0.17252942, 0.17164107, 0.16064716, 0.14564287, 0.14578649,
       0.14152733, 0.1354919 , 0.11175379, 0.1380746 , 0.12547517,
       0.15136653, 0.13984282, 0.18308302, 0.12271885, 0.15289988,
       0.13492309, 0.13499516, 0.13373476, 0.1034279 , 0.14278288,
       0.14574681, 0.11614764, 0.11256923, 0.14796558, 0.11459825,
       0.12417535, 0.15693744, 0.14159134, 0.11885544, 0.13164357,
       0.13445257, 0.13527885, 0.13472062, 0.12027512, 0.12072214,
       0.15361264, 0.12973932, 0.11003032, 0.13575847, 0.11980422,
       0.1187932 , 0.11152574, 0.14656588, 0.13885414, 0.13960315,
       0.12921241, 0.09522926, 0.14543513, 0.14980696, 0.11318417,
       0.10785905, 0.13858491, 0.11922434, 0.11760534, 0.12059705,
       0.12150726, 0.1184712 , 0.11084933, 0.10894509, 0.10107464,
       0.10258616, 0.1094653 , 0.095096  , 0.10059849, 0.10931144,
       0.11704954, 0.12639652, 0.13283708, 0.10203757, 0.10787873])

def l(para,args):
    xi,sigma = para
    y = args
    k = len(y)
    SUM = log(1+(xi*y)/sigma)
    return (-k*log(sigma)-(1+xi)/xi*np.sum(SUM))

def constrain1(x):
    xi,sigma = x
    return sigma
def constrain2(x):
    xi,sigma = x
    return -xi
def constrain3(x,*args):
    xi,sigma = x
    yi = args
    return 1+xi*yi/sigma
    

con1 = {'type': 'ineq', 'fun': constrain1}
con2 = {'type': 'ineq', 'fun': constrain2}
con3 = {'type': 'ineq', 'fun': constrain3,'args': y}
cons = [con1,con2,con3]

x0 =(-0.1,2)   
bound =((-np.inf,None),(None,np.inf))
res = minimize(l,x0,(y),bounds = bound,constraints=cons)

我遇到的第一个问题是实现约束 3，即对数内的项对于每个数据点都大于零 (1+(xi y_i)/ sigma>0）。我试图在 for 循环中添加约束

con1 = {'type': 'ineq', 'fun': constrain1}
con2 = {'type': 'ineq', 'fun': constrain2}
cons = [con1,con2]
for yi in y:
    con3 = {'type': 'ineq', 'fun': constrain3,'args': yi}
    cons.append(con3)

但这也没有用。

我的第二个问题是约束只支持 >= 或 =。有没有办法只实现 >，尤其是不相等 (!=)？我想要 xi != 0 和 sigma > 0。我很惊讶我没有找到任何关于该主题的内容，因为这应该是数学中的标准问题，或者我看起来不够努力

我在 mathlab 中找到了一种解决方法，我使用 xi 进行了一项优化，另一项使用 xi>=0 并选择了更合适的。 sigma 似乎并没有受到它可能为零的约束的太大影响。但我更喜欢不需要在不同解决方案之间进行评估的解决方案。

【问题讨论】：

一些评论：严格的不等式在理论上大多数优化框架中不受支持并且没有多大意义。需要显式的基于 epsilon 的公式 (>= 0 + eps)，并且您必须确定 epsilon 是 10^-100 还是 10^-2。 scipy.minimize 主要是关于数值优化（在最佳情况下：凸），但 x != 0 将空间划分为非凸集/析取，如果没有分支/离散优化技术，您将无法做到这一点。简而言之：您要求的东西可以让我们轻松解决所有这些 NP 难题。必须有一个陷阱。
选择的路径通常取决于用例，一些知识肯定会有所帮助。如果没有上下文并且仅通过查看您的方程式，我只想说，x != 0 is 约束是一种避免被零除的技巧，并且永远不应作为分支不等式处理。在 scipy 优化中，您可能会抓住这种情况并做一些类似次梯度的事情。不过，如果这行得通，那就是另一回事了。
感谢@sascha 的解释。你说的对。听起来很简单，不会有问题。我尝试在统计极值评估的背景下将通用帕累托分布拟合到我的数据中，并且约束 `xi!=0` 确实存在，因此您不能除以零。这是我第一次尝试使用 scipy.minimize ，我的数学背景仅限于工程观点。我应该如何在我的代码中实现它？
我会忽略它并且不使用初始启动，其中该值 not 设置为零。它变为零的机会非常渺茫（甚至可以在你的函数中处理一些 if；在那里失去平滑度希望不会杀死求解器）。您必须进行一些实验。我可以想象你会遇到不止一个数值问题，包括需要对 sigma 进行 epsilon 调整。此外，如果它是非凸的（不确定），您可能也想使用多头。

标签： python optimization scipy statistics minimize

【解决方案1】：

import numpy as np

from numpy import log
from scipy.optimize import minimize

y = np.array([0.21863326, 0.19805154, 0.22953017, 0.21906749, 0.22067327,
              0.20638931, 0.17845443, 0.20008429, 0.21702199, 0.16334912,
              0.18480577, 0.17172182, 0.16495525, 0.15907978, 0.17029919,
              0.14020628, 0.16528562, 0.17141436, 0.14978351, 0.1329871,
              0.14036109, 0.15894933, 0.16783223, 0.17372222, 0.15986161,
              0.1654368, 0.16348146, 0.15595923, 0.15192792, 0.12272897,
              0.17252942, 0.17164107, 0.16064716, 0.14564287, 0.14578649,
              0.14152733, 0.1354919, 0.11175379, 0.1380746, 0.12547517,
              0.15136653, 0.13984282, 0.18308302, 0.12271885, 0.15289988,
              0.13492309, 0.13499516, 0.13373476, 0.1034279, 0.14278288,
              0.14574681, 0.11614764, 0.11256923, 0.14796558, 0.11459825,
              0.12417535, 0.15693744, 0.14159134, 0.11885544, 0.13164357,
              0.13445257, 0.13527885, 0.13472062, 0.12027512, 0.12072214,
              0.15361264, 0.12973932, 0.11003032, 0.13575847, 0.11980422,
              0.1187932, 0.11152574, 0.14656588, 0.13885414, 0.13960315,
              0.12921241, 0.09522926, 0.14543513, 0.14980696, 0.11318417,
              0.10785905, 0.13858491, 0.11922434, 0.11760534, 0.12059705,
              0.12150726, 0.1184712, 0.11084933, 0.10894509, 0.10107464,
              0.10258616, 0.1094653, 0.095096, 0.10059849, 0.10931144,
              0.11704954, 0.12639652, 0.13283708, 0.10203757, 0.10787873])


def l(para, args):
    xi, sigma = para
    y = args
    k = len(y)
    # Define a list of functions to use np.sum() on
    SUM = [log(1 + (xi * y_i) / sigma) for y_i in y]
    return -k * log(sigma) - (1 + xi) / xi * np.sum(SUM)


def constrain1(x):
    (xi, sigma) = x
    return sigma


def constrain2(x):
    xi, sigma = x
    return -xi


# Use a class to prevent the python garbage collector to kick in
class Constraint3:
    def __init__(self, y_i):
        self.y_i = y_i

    def make_constraint(self, x):
        xi, sigma = x
        return 1 + xi * self.y_i / sigma


con1 = {'type': 'ineq', 'fun': constrain1}
con2 = {'type': 'ineq', 'fun': constrain2}
constraints_list = [con1, con2]

for y_i in y:
    # Make each constraint function unique by instantiating the above class with the
    # new y_i value. The usage of individually defined function is not advised, because
    # some of them could get garbage collected, when there are too many items in the
    # array.
    constraints_list.append({'type': 'ineq', 'fun': Constraint3(y_i).make_constraint})

# Use proper ndarray-type
x0 = np.array([-0.1, 2])
bound = ((-np.inf, None), (None, np.inf))

res = minimize(l, x0, y, bounds=bound, constraints=constraints_list,
               options={'disp': True})
print("            Parameter values:")
print("              xi:   ", res.x[0])
print("              sigma:", res.x[1])

# Optimization terminated successfully    (Exit mode 0)
#         Current function value: -1739.2412627695364
#         Iterations: 24
#         Function evaluations: 82
#         Gradient evaluations: 24
#         Parameter values:
#           xi:    -1928.571453953235
#           sigma: 35762853.400437616

【讨论】：