Nichtlineare Anpassung

Jens Wagner 07/2015

Download this notebook
Download dispersion.txt
Informationen zum Fitten

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import pylab as py 
%matplotlib inline
In [2]:
lamb, winkel, err = py.loadtxt('data/dispersion.txt', unpack=True, skiprows=1) #Messwerte importieren
# Versuch Prismenspektrometer
#lamb: Wellenlaenge lambda, winkel: Ablenkwinkek, err: Fehler des Winkels
plt.figure(figsize=(8,6))
plt.errorbar(lamb, winkel, yerr=err, fmt='.') 
plt.xlabel('Wellenlaenge / nm')
plt.ylabel('Winkel / deg')
plt.title('Dispersion')

plt.show()
In [3]:
#Fitfunktion360/np.pi*np.arcsin(0.5*(n0+k/(lamb-l))) - 60
def fitFunc(lamb, n0, k, l):
    return 360/np.pi*np.arcsin(0.5*(n0+k/(lamb-l))) - 60
In [4]:
init_vals = [1.8, 8, 200] #Fitparameter initialisieren (n0,k,l)
#Falls der Fit nicht konvergiert muessen die Startparameter geaendert werden
fitParams, fitCovariances = curve_fit(fitFunc, lamb, winkel ,p0=init_vals, sigma=err)
#curve_fit() uebergibt den gefundenen Parametersatz und die Kovarianzmatrix
In [5]:
print(fitParams)
print(fitCovariances) #Fehler der Parameter sind die Wurzeln der Diagonalelemente der Kovarianzmatrix

[   1.80535046    8.49799329  208.58402713]
[[  1.17439415e-06  -6.73178942e-04   1.04805525e-02]
 [ -6.73178942e-04   3.92209488e-01  -6.17979831e+00]
 [  1.04805525e-02  -6.17979831e+00   9.84944832e+01]]
In [6]:
print("n0=", fitParams[0], ", Standardfehler=", np.sqrt(fitCovariances[0][0]))
print("k=" , fitParams[1], ", Standardfehler=", np.sqrt(fitCovariances[1][1]))
print("l=" , fitParams[2], ", Standardfehler=", np.sqrt(fitCovariances[2][2]))

n0= 1.80535045755 , Standardfehler= 0.00108369467435
k= 8.49799329294 , Standardfehler= 0.626266307641
l= 208.584027134 , Standardfehler= 9.92443868534
In [7]:
#Diagramm mit angepasster Kurve zeichnen
x=np.linspace(400,700,50) #x-Werte fuer die Fitfunktion 
plt.figure(figsize=(8,6))
plt.plot(x, fitFunc(x, fitParams[0], fitParams[1], fitParams[2]), label=r'$360/\pi\, \arcsin(1/2(n_0+k/(\lambda-l))) - 60$')
plt.errorbar(lamb, winkel, yerr=err, fmt='.', color='b', label=r'Messwerte') 
plt.xlabel('Wellenlaenge / nm')
plt.ylabel('Winkel / deg')
plt.legend(loc='best',prop={'size':16})
plt.title('Dispersion')
Out[7]:
<matplotlib.text.Text at 0x77cfa20>
In [8]:
#Chi_quadrat und Chi_quadrat_reduziert berechnen
chisquare=np.sum((fitFunc(lamb, fitParams[0], fitParams[1], fitParams[2])-winkel)**2/err**2)
dof=len(winkel)-len(init_vals) #dof: degees of freedom, Zahl der Messwerte abzueglich Zahl der Parameter
print("dof=",dof)
print("chisquare=", chisquare, ", chisquare_red=",chisquare/dof)

dof= 7
chisquare= 7.6898070518 , chisquare_red= 1.09854386454

Goodness of fit

Berechne die Wahrscheinlichkeit einen chi2 Wert zu bekommen, der groesser oder gleich dem eben berechneten chi2 Wert ist

In [9]:
from scipy.stats import chi2
print("propability="+str(round(100-chi2.cdf(chisquare,dof)*100,2))+"%") #cdf: cumulative distribution function

propability=36.07%
In [10]:
#Residuen in Einheiten des Standardfehlers
residuen=(winkel-fitFunc(lamb, fitParams[0], fitParams[1], fitParams[2]))/err
plt.plot(lamb, residuen, marker='o', linewidth=0)
Out[10]:
[<matplotlib.lines.Line2D at 0x79fec88>]
In [11]:
#plt.hist(residuen) #zu wenig Daten! nicht sehr aussagekraeftig