最小化使用数据拟合蟒函数

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

Unconcsiosely，这个问题重复我的其他问题。正确的答案是：

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) 

Answer 2

这里是我的问题的理解。 首先，我通过下面的代码

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协方差矩阵是简化形式，需要重新调整比例。