使用Python最小化函数进行数据拟合

Question

我有一个功能如下

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

是积分范围。 我需要您的帮助，感谢您的时间和关注。 请不要评价负值。我需要你的答案。

Answer 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

，协方差矩阵采用简化形式，需要重新缩放。

Answer 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)