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

Question

我尝试通过最似然法将统计方程的参数获取到我的数据中。我想适应约束的方程看起来像这样： 方程和约束

最佳拟合是在我的似然函数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<=0进行了一次优化，使用xi>=0进行了另一次优化，并选择了更好的拟合。sigma似乎并没有受到它可能为零的约束的太大影响。但我更喜欢一个解决方案，我不必在不同的解决方案之间进行评估。

标签： pythonoptimizationscipystatisticsminimize

Answer 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