# 数据分析--拟合

1.求根

newton法

F(X)的雅克比矩阵

2.最小二乘拟合

In 1929 Edwin Hubble examined the data in Table relating the radial velocity v of 24 extra galactic nebulae to their distance r from our galaxy. Although there was considerable scatter in the data, he fit them with a straight line:

y(x) = a0+a1 x
a0 = 0.032 a1 = 0.9

Measured total cross sections in millibarns of neutron-carbon scattering as a function of neutron energy in KeVs.

import math
import numpy as np
from numpy import *
import matplotlib.pyplot as plt
import matplotlib

plt.figure(figsize=(20,16), dpi=100)
plt.rc('font',family='Times New Roman')
matplotlib.rcParams['font.family']='STSong'
plt.rcParams['axes.unicode_minus']=False
matplotlib.rcParams['font.size']=20

ND=9
num =3
E=np.array([0, 13 ,25 ,38, 50, 63, 75, 87, 100])
S=np.array([10.6, 16.0, 45.0, 83.5, 52.8, 19.9, 10.8, 8.25, 4.7])
sig = np.ones_like(S)
#sig =np.array([9.34, 17.9, 41.5, 85.5, 51.5, 21.5, 10.8, 6.29, 4.14])    # Error bars
x = np.array([20000,40,220],dtype=float)
plt.plot(E, S, 'bo',markersize=20,alpha=0.4,label="采样数据点")

plt.title('Least-Squares Fit of linear to Blue Data',fontsize=20)
plt.xlabel('E',fontsize=32);  plt.ylabel('$\sigma$',fontsize=32);   plt.grid(True)
def Fun(x,num):
F = np.zeros((num),dtype=float)
for j in range(ND):
F[0] += (S[j]-x[0]/((E[j]-x[1])**2+x[2]))/sig[j]**2/((E[j]-x[1])**2+x[2])
F[1] += -2*x[0]*(E[j]-x[1])*(S[j]-x[0]/((E[j]-x[1])**2+x[2]))/sig[j]**2/((E[j]-x[1])**2+x[2])**2
F[2] += -x[0]*(S[j]-x[0]/((E[j]-x[1])**2+x[2]))/sig[j]**2/((E[j]-x[1])**2+x[2])**2

return F

def dfun(x,num):
df = np.zeros((num,num),dtype=float)
dx = 0.00001                           #
x1_plus = np.copy(x)
x1_minus = np.copy(x)

for i in range(0,num):
for j in range(0,num):
x1_plus = np.copy(x)
x1_minus = np.copy(x)
x1_plus[j] = x1_plus[j]+dx/2           #x+dx
x1_minus[j] = x1_minus[j]-dx/2
df[i,j] = (Fun(x1_plus,num)[i]-Fun(x1_minus,num)[i])/dx   #f(x+dx)-f（x）/dx
df_1 = np.linalg.inv(df)                              #计算逆矩阵
return df_1

def Newton(x,num):
x1 = np.copy(x)
i = 0
delta = np.copy(x)
#    dfun0=dfun(x,num)
while(np.sum(abs(delta)) > 1.e-6 and i < 100):  #控制循环次数
x1 = x-np.dot(dfun(x,num),Fun(x,num))  #公式
delta = x1-x
x = x1
i = i+1
print(x,)
return x

a = Newton(x,num)
print(a)

EX=np.linspace(0,100,1000)
G_E=x[0]/((EX-x[1])**2+x[2])
G = a[0]/((EX-a[1])**2+a[2])
points = a[0]/((E-a[1])**2+a[2])
plt.plot(EX, G_E,'k', linewidth=3.0,label="初猜曲线")
plt.plot(EX, G,'r:',linewidth=5.0,label="最终拟合曲线")
plt.plot(E, points, 'ro',markersize= 20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.legend()
plt.show()

results

[ 14956.38460591 38.99233378 177.84417855]

[ 16568.5542651 38.80318039 196.28768571]

[ 16786.04125645 38.81295456 198.89084396]

[ 16789.85161833 38.81293923 198.94072711]

[ 16789.85282055 38.81293926 198.94074454]

[ 16789.85282055 38.81293926 198.94074454]

[ 16789.85282055 38.81293926 198.94074454]

import numpy as np
from scipy.optimize import leastsq
import pylab as pl

def func(x, p):
a1, a2, a3 = p
return a1/((x-a2)**2+a3)

def residuals(p, y, x):
return y - func(x, p)

x = np.linspace(0, 100, 1000)

p0 = [16000,40, 200, ]
x0=np.array([0, 13 ,25 ,38, 50, 63, 75, 87, 100])
y0=np.array([10.6, 16.0, 45.0, 83.5, 52.8, 19.9, 10.8, 8.25, 4.7])
y1 = y0 + 2 * np.random.randn(len(x0))

plsq = leastsq(residuals, p0, args=(y1, x0))

print (u"拟合参数", plsq[0])
plt.title(u'洛伦兹函数最小二乘法拟合',fontsize=20)
pl.plot(x0, y0, 'bo',markersize=20,alpha=0.4, label=u"真实数据")
pl.plot(x0, y1, 'y*', markersize=20,alpha=0.4, label=u"带噪声的实验数据")
pl.plot(x, func(x, plsq[0]), 'r:',linewidth=5.0, label=u"拟合曲线")
plt.plot(x, func(x, p0),'k', linewidth=3.0,label="初猜曲线")
plt.plot(x0, func(x0, plsq[0]), 'ro',markersize= 20)
plt.xlabel('E',fontsize=32);  plt.ylabel('$\sigma$',fontsize=32);   plt.grid(True)
pl.legend()
pl.show()

https://blog.sciencenet.cn/blog-566204-1173092.html

## 全部精选博文导读

GMT+8, 2024-9-7 20:25