python - 最小化非线性方程组（指数积分）

Question

一般的：

我正在使用最大熵来查找正整数向量的分布，我可以估计均值和方差，并尝试找到三个方程 a 和 b，

方程：

1. 积分(exp(a*x^2+bx+c) from (0 , infinity))-1

2. 积分(x exp(a x^2+bx+c)from (0 , infinity))- 均值

3. 积分(x^2*exp(a*x^2+bx+c) from (0 , infinity))- mean^2 - var

（[0,∞) 之间的积分）

问题：

我正在尝试使用数值求解器，并且我使用了 sympy 的 fsolve 但我想我缺少一些知识。

我的代码：

import numpy as np
import sympy as sym
from scipy.optimize import *


def myFunction(x,*data):
    y = sym.symbols('y')
    m,v=data
    F = [0]*3
    x[0] = - abs(x[0])
    print(x)
    F[0] = (sym.integrate(sym.exp(x[0] * y ** 2 + x[1] * y + x[2]), (y, 0,sym.oo)) -1).evalf() 
    F[1] = (sym.integrate(y*sym.exp(x[0] * y ** 2 + x[1] * y + x[2]), (y, 0,sym.oo))-m).evalf()
    F[2] = (sym.integrate((y**2)*sym.exp(x[0] * y ** 2 + x[1] * y + x[2]), (y,0,sym.oo)) -v-m).evalf() 
    print(F)
    return F


data = (10,3.5) # mean and var for example
xGuess = [1, 1, 1]
z = fsolve(myFunction,xGuess,args = data)
print(z)

我的结果不是那么准确，有没有更好的方法来解决它？

1. 积分(exp(a*x^2+bx+c))-1 = 5.67659292676884

2. 积分(x exp(a x^2+bx+c))- 均值 = −1.32123173796713

3. 积分(x^2*exp(a*x^2+bx+c))- 平均^2 - var = −2.20825624606312

谢谢

Answer 1

我已经sympy用numpy和 lambdas（内联函数）重写了问题。另请注意，在您的问题陈述中，您用 $mean^2$ 减去第三个等式，但在您的代码中，您只减去 $mean$。

import numpy as np
from scipy.optimize import minimize
from scipy.integrate import quad


def myFunction(x,data):
    m,v=data
    F = np.zeros(3)  # use numpy array

    # use scipy.integrade.quad for integration of lambda functions
    # quad output is (result, error), so we just select the result value at the end
    F[0] = quad(lambda y: np.exp(x[0] * y ** 2 + x[1] * y + x[2]), 0, np.inf)[0] -1
    F[1] = quad(lambda y: y*np.exp(x[0] * y ** 2 + x[1] * y + x[2]), 0, np.inf)[0] -m
    F[2] = quad(lambda y: (y**2)*np.exp(x[0] * y ** 2 + x[1] * y + x[2]), 0, np.inf)[0] -v-m**2

    # minimize the squared error
    return np.sum(F**2)


data = (10,3.5) # mean and var for example
xGuess = [-1, 1, 1]
z = minimize(lambda x: myFunction(x, data), x0=xGuess,
             bounds=((None, 0), (None, None), (None, None)))  # use bounds for negative first coefficient
print(z)

# x: array([-0.99899311,  2.18819689,  1.85313181])

这似乎更合理？