scipy中的数学范围错误最小化

Question

我想用最大似然法将六个参数拟合到非常丑陋的分布函数中。为此，我尝试使用scipy.optimize.minimize。

这是一段代码

import math
form scipy.optimize import minimize
import numpy as np

#generate some data
xdata = np.random.lognormal(0,1,812)
#function for the Log likelihood
def mfpdf2(params):
    c1 = params[0]
    A1 = params[1]
    a1 = params[2]
    c2 = params[3]
    A2 = params[4]
    a2 = params[5]

    LL_vec = [math.log(c1*(math.exp(-A1*x) - math.exp(-a1*x))+c2*(math.exp(-A2*x) - math.exp(-a2*x))) for x in xdata]
    LL = -sum(LL_vec)


return LL

#try to find max likelihood (minimize negative loglikelihood)

start_params = [1,1,2,1,1,2]

pars = minimize(mfpdf2, start_params)

这段代码返回一个错误：

File "C:/Users/Robert/Desktop/python/pokusy_analyza_multiexp.py", line 79, in <listcomp>
LL_vec = [math.log(c1*(math.exp(-A1*x) - math.exp(-a1*x))+c2*(math.exp(-A2*x) - math.exp(-a2*x))) for x in xdata]

ValueError: math domain error

我做错了什么？

Answer 1

看起来你可能正在记录一个未定义的负数

Answer 2

如果x 为负数，您在math.log 函数中的表达式将返回一个负数。负数的对数只为复数定义，所以常规的对数函数会给你一个ValueError: math domain error。您可以将x_data 转换为复数并使用np.log，如果这将在您的特定应用程序中从minimize 产生可接受的结果：

xdata = np.random.lognormal(0,1,812).astype(np.complex)
...
LL_vec = [np.log(c1*(math.exp(-A1*x) - m ...

或者您可以指定不同的最小化技术（并非所有都支持定义的界限）并指定值的界限。这不支持像 A1 > a1 这样的符号边界，因此您必须重新排序变量以建立关系：

from scipy.optimize import minimize
import numpy as np

xdata = np.random.lognormal(0,1,812)

def mfpdf2(params):

    c1 = params[0]
    A1 = params[1]
    offset1 = params[2] #this can obviously be condensed to one line
    a1 = A1 + offset1 #bound offset to be positive so a is always > A 
    c2 = params[3]
    A2 = params[4]
    offset2 = params[5]
    a2 = A2 + offset2

    LL_vec = [np.log(
                     c1*(np.exp(-A1*x) - np.exp(-a1*x))+
                     c2*(np.exp(-A2*x) - np.exp(-a2*x))
                 ) for x in xdata]
    #this will not account for the possibility that 
    #   one of c1*(...) or c2*(...) is negative but 
    #   the sum is still positive. This could conceiveably 
    #   be achieved with more ratios or offsets instead of 
    #   direct values, but would make the math real nasty.
    LL = -sum(LL_vec)
    return LL

start_params = [1,1,2,1,1,2]
pars = minimize(mfpdf2, 
                start_params, 
                method='L-BFGS-B',
                #now define your bounds
                bounds=((None, None), # c1
                        (None, None), # A1
                        (0, None),    # a1 - A1 (makes a1 always larger than A1)
                        (None, None), # c2
                        (None, None), # A2
                        (0, None)))   # a2