使用python最小化函数进行数据拟合答案

【问题标题】：Minimizing a function using python for data fitting使用python最小化函数进行数据拟合
【发布时间】：2017-08-14 14:58:09
【问题描述】：

我有如下功能

q = 1 / sqrt( ((1+z)**2 * (1+0.01*o_m*z) - z*(2+z)*(1-o_m)) )
h = 5 * log10( (1+z)*q ) + 43.1601

我有上述方程的实验答案，一旦我必须将一些数据放入上述函数并求解下面的方程

chi=(q_exp-q_theo)**2/err**2  # this function is a sigma, sigma chi from z=0 to z=1.4 (in the data file)

z、err 和 q_exp 在数据文件（2.txt）中。现在我必须为o_m (0.2 to 0.4) 选择一个范围，然后找到o_m 的范围，chi 函数将被最小化。

我的代码是：

from math import *
from scipy.integrate import quad

min = None
l = None
a = None
b = None
c = 0


def ant(z,om,od):
    return 1/sqrt( (1+z)**2 * (1+0.01*o_m*z) - z*(2+z)*o_d )


for o_m in range(20,40,1):
    o_d=1-0.01*o_m
    with open('2.txt') as fp:
        for line in fp:
            n = list( map(float, line.split()) )
            q = quad(ant,n[0],n[1],args=(o_m,o_d))[0]
            h = 5.0 * log10( (1+n[1])*q ) + 43.1601
            chi = (n[2]-h)**2 / n[3]**2
            c = c + chi
        if min is None or min>c:
            min = c
            l = o_m                            
        print('chi=',q,'o_m=',0.01*l)

n[1],n[2],n[3],n[4]分别在数据文件中是z1,z2,q_exp和err。 z1 和 z2 是积分范围。我需要您的帮助，感谢您的时间和关注。请不要评价负值。我需要你的答案。

【问题讨论】：

1：你的问题是什么？ 2：请分享一些最小的数据集。 3：为什么ant()有o_m和o_d，而上面的q只有o_m。
为什么不使用scipy.optimize.leastsq？ BTW，如果数据不是太大，一开始就加载一次，可能用numpy.loadtxt()
此外，还有一些错别字；在最后一个print() 中，您使用' 打开但使用" 关闭，您可能想要打印c 而不是q，您应该命名，例如c2 因为它实际上是正方形。这样可以避免混淆。有些缩进似乎是错误的。评论输入 pythonic：min == None 有效，但 min is None 看起来更好。甚至可能是if not min。您真的想与h 或c 进行比较吗？
这是卡方检验。我们有一个集成，其中将使用数据文件来获得最终答案。我们必须最小化卡方并找到使卡方最小化的参数。 om 和 od 是 print 中 main 函数中的 omega d 和 omega m 我们只需要 om 或 od。
我更正了你提到的一些部分

标签： python data-fitting integrate quad

【解决方案1】：

这是我对问题的理解。首先我通过以下代码生成一些数据

import numpy as np
from scipy.integrate import quad
from random import random


def boxmuller(x0,sigma):
    u1=random()
    u2=random()
    ll=np.sqrt(-2*np.log(u1))
    z0=ll*np.cos(2*np.pi*u2)
    z1=ll*np.cos(2*np.pi*u2)
    return sigma*z0+x0, sigma*z1+x0


def q_func(z, oM, oD):
    den= np.sqrt( (1.0 + z)**2 * (1+0.01 * oM * z) - z * (2+z) * (1-oD) )
    return 1.0/den


def h_func(z,q): 
    out = 5 * np.log10( (1.0 + z) * q ) + .25#43.1601
    return out


def q_Int(z1,z2,oM,oD):
    out=quad(q_func, z1,z2,args=(oM,oD))
    return out

ooMM=0.3
ooDD=1.0-ooMM

dataList=[]
for z in np.linspace(.3,20,60):
    z1=.1+.1*z*.01*z**2
    z2=z1+3.0+.08+z**2
    q=q_Int(z1,z2,ooMM,ooDD)[0]
    h=h_func(z,q)
    sigma=np.fabs(.01*h)
    h=boxmuller(h,sigma)[0]
    dataList+=[[z,z1,z2,h,sigma]]
dataList=np.array(dataList)

np.savetxt("data.txt",dataList)

然后我将按照以下方式进行调整

import matplotlib
matplotlib.use('Qt5Agg')

from matplotlib import pyplot as plt
import numpy as np
from scipy.integrate import quad
from scipy.optimize import leastsq 

def q_func(z, oM, oD):
    den= np.sqrt( (1.0 + z)**2 * (1+0.01 * oM * z) - z * (2+z) * (1-oD) )
    return 1.0/den


def h_func(z,q): 
    out = 5 * np.log10( (1.0 + z) * q ) + .25#43.1601
    return out


def q_Int(z1,z2,oM,oD):
    out=quad(q_func, z1,z2,args=(oM,oD))
    return out


def residuals(parameters,data):
    om,od=parameters
    zList=data[:,0]
    yList=data[:,3]
    errList=data[:,4]
    qList=np.fromiter( (q_Int(z1,z2, om,od)[0] for  z1,z2 in data[ :,[1,2] ]), np.float)
    hList=np.fromiter( (h_func(z,q) for z,q in zip(zList,qList)), np.float)
    diffList=np.fromiter( ( (y-h)/e for y,h,e in zip(yList,hList,errList) ), np.float)
    return diffList

dataList=np.loadtxt("data.txt")

###fitting 
startGuess=[.4,.8]
bestFitValues, cov,info,mesg, ier = leastsq(residuals, startGuess , args=( dataList,),full_output=1)
print bestFitValues,cov

fig=plt.figure()
ax=fig.add_subplot(1,1,1)
ax.plot(dataList[:,0],dataList[:,3],marker='x')


###fitresult
fqList=[q_Int(z1,z2,bestFitValues[0], bestFitValues[1])[0] for z1,z2 in zip(dataList[:,1],dataList[:,2])]
fhList=[h_func(z,q) for z,q in zip(dataList[:,0],fqList)]
ax.plot(dataList[:,0],fhList,marker='+')


plt.show()

输出

>>[ 0.31703574  0.69572673] 
>>[[  1.38135263e-03  -2.06088258e-04]
>> [ -2.06088258e-04   7.33485166e-05]]

和图表请注意，对于leastsq，协方差矩阵是简化形式，需要重新调整。

【讨论】：

非常感谢，非常完美。但拟合数据不是我所期望的。我们更改 om 并检查输出是否存在 om 使 chi 函数最小化。我的代码是正确的。但需要一些修改
@Ehsan 我不会称之为“正确”的原始版本；无论如何，这就是为什么我在上面问你真正想要的是什么。您以低效的方式显示一维拟合。正如您拥有o_m 和o_d 一样，我介绍了一般的2D 案例，您可以轻松地将其修改为1D。如果你想要一个可以工作的版本，你真的应该发布一个最小的数据集，也许是 20 到 30 个点。除此之外，我的代码正在执行您的代码正在执行的操作，除了您正在使用的低效 for 循环。从这个角度来看，它是您的代码的修订版，就像修订版 10 省略了步骤 2 到 9。
@Ehsan 我真的建议真正回答我在上面的 cmets 中提出的问题。
非常感谢。我从###fitting 到最后一遍又一遍地阅读它，但我不太明白它们是什么。以及它们与我的问题的最后一部分之间的关系是什么。事实上，我再次修改了我的代码，但问题是：我必须找到一个导致最小“chi”或最小“c”的“o_m”。我可以找到它，但方法错误。在我在 for 循环中为 o_m 选择的每个范围内，该循环的最小值导致“c”的最小值。例如，在我的代码中，o_m 有 (20,40,1)，最少有 20 个结果。这是错误的，我不知道如何解决它
因为最小值必须介于 25 到 31 之间，或者不是最小数字或最大数字或我们在循环中的任意范围的值。

【解决方案2】：

不自觉地，这个问题与我的另一个问题重叠。正确答案是：

 from math import *
 import numpy as np
 from scipy.integrate import quad
 min=l=a=b=chi=None
 c=0
 z,mo,err=np.genfromtxt('Union2.1_z_dm_err.txt',unpack=True)
 def ant(z,o_m):            #0.01*o_m  is steps of o_m
     return 1/sqrt(((1+z)**2*(1+0.01*o_m*z)-z*(2+z)*(1-0.01*o_m)))
 for o_m in range(20,40):
     c=0
     for i in range(len(z)):
         q=quad(ant,0,z[i],args=(o_m,))[0]     #Integration o to z
         h=5*log10((1+z[i])*(299000/70)*q)+25     #function of dL
         chi=(mo[i]-h)**2/err[i]**2               #chi^2 test function
        c=c+chi
        l=o_m
        print('chi^2=',c,'Om=',0.01*l,'OD=',1-0.01*l)

【讨论】：