1 | import os |
2 | import numpy as np |
3 | import matplotlib.pyplot as plt |
4 | from scipy.optimize import curve_fit |
5 | os.chdir(r'D:\data\scripts\wordpress\ex55') |
6 | os.getcwd() |
7 | x0,A,gamma=12,3,5 |
8 | n=200 |
9 | x=np.linspace(1,20,n) |
10 | yexact=A*gamma**2/(gamma**2+(x-x0)**2) |
11 | #Learning Scientific Programming with Python |
12 | # Add some noise with a sigma of 0.5 apart from a particularly noisy region |
13 | # near x0 where sigma is 3 |
14 | sigma=np.ones(n)*0.5 |
15 | sigma[np.abs(x-x0+1)<1]=3 |
16 | noise=np.random.randn(n)*sigma |
17 | y=yexact+noise |
18 |
19 | def f(x,x0,A,gamma): |
20 | #The Lorentzian entered at x0 with amplitude A and HWHM gamma |
21 | return A*gamma**2/(gamma**2+(x-x0)**2) |
22 |
23 | def rms(y,yfit): |
24 | return np.sqrt(np.sum((y-yfit)**2)) |
25 |
26 | # Unweighted fit |
27 | p0=10,4,2 |
28 | popt,pcov=curve_fit(f,x,y,p0) |
29 | yfit=f(x,*popt) |
30 | print('Unweighted fit parameters:',popt) |
31 | print('Covarian cematrix:') |
32 | print(pcov) |
33 | print('rms error in fit:',rms(yexact,yfit)) |
34 | print() |
35 | # Weighted fit |
36 | popt2,pcov2=curve_fit(f,x,y,p0,sigma=sigma,absolute_sigma=True) |
37 | yfit2=f(x,*popt2) |
38 | print('Weighted fit parameters:',popt2) |
39 | print('Covariancematrix:') |
40 | print(pcov2) |
41 | print('rms error in fit:',rms(yexact,yfit2)) |
42 | plt.plot(x,yexact,label='Exact') |
43 | plt.plot(x,y,'o',label='Noisy data') |
44 | plt.plot(x,yfit,label='Unweighted fit') |
45 | plt.plot(x,yfit2,label='Weighted fit') |
46 | plt.ylim(-1,4) |
47 | plt.legend(loc='lower center') |
48 |
49 | plt.tight_layout() |
50 | plt.savefig("example55.png", dpi=100) |
51 | plt.show() |
52 | plt.close() |
Recent Comments