Embedded Files
Skal du bruge koden til eksponentiel regression så scroll ned til "Her er koden til exponentiel regression"
Exponentiel regression
Exponentiel regression
Her skal vi kikke på fit som bruger nedenstående funktion.
Til højre ser du et billede af en
y = a * exp(-b*x) + c
eksponentiel regression. Lige som lineær regression (fx i excel hvor du tilføjer en tendenslinje) prøver eksponentil regression at finde den bedste eksponentielle funktion igennem de givne datapunkter.
Om at definere en funktion
Om at definere en funktion
Prøv at køre koden til højre og se hvad der sker. Funktionens navn volumen har du selv valgt og den "returnerer" det som du angiver den skal gøre.
def volumen(laengde,bredte,hoejde):
return laengde*bredte*hoejde
print(volumen(2,3,4))
Tjek nu koden længere nede. Der defineres en funktion func som beregner y-værdierne til en e-funktion som tager en masse x-værdier som input. a,b,c er konstanter som definerer funktionen.
Funktionen bliver "kaldt" til routinen curve_fit som har x_vaerdier, y_vaerdier som input og p0=[250, 0.13, 0] som startværdier på regressionsberegningen.
def func(x, a, b, c):
return a * np.exp(-b * x) + c
popt, pcov = curve_fit(func, x_vaerdier, y_vaerdier, p0=[250, 0.13, 0])
Til en lineær model skulle funktionen se således ud
def func(x, a, b):
return a * x + b
popt, pcov = curve_fit(func, x_vaerdier, y_vaerdier, p0=[1, 0])
Kommentar
Kommentar
scipy.optimize er den pakke som indeholder rutinen "curve_fit" som bruges her igen lige som i det lineære fit.
x_vaerdier og y_vaerdier har jeg selv fundet på så jeg har noget at fitte på. De skal selvfølgelig erstattes af "rigtige" værdier.
Denne her rutine bliver brugt i radioaktivitetsprojektet til at finde halveringstiden på udvalgte henfald.
Jeg har skrevet nogle kommentarer med rød ind i programmet
Som i kender det måske fra LoggerPro skal fitteprogrammer have "startværdier" til at prøve at finde den rigtige funktion. Det skal curve_fit også. Startværdierne er
p0=[250, 0.5, 0])
Dem har jeg gættet mig til. Hvordan?
Funktionen er y = a*exp(-b x) + c
Startværdier for a b og c
Startværdier for a b og c
p0=[0] = 250 fordi den er tæt på det støreste værdi i mit datasæt
p0=[1] = 0.13 fordi det er min "henfaldskonstant" Vel ... et hurtigt overslag giver at halveringstiden fra 250 til 120 er ca. 5s (x-værdierne) ln(2)/5 = 0,13 som er et godt bud på en b-værdi.
p0=[2] = 0 c-værdien regner jeg med er 0 fordi funktionen starter i 251.
Python er ikke dum og jeg har prøvet lidt med forskellige værdier og den er VIRKELIG large med at finde et fit. Den er klart mest følsom overfor forkerte b-værdier.
Se hvor fin man nu kan beregne halveringtiden.
Som antydet ovenfor er det IKKE trivielt at skulle finde startværdieerne for fittet, her skal man lige kikke på de fysiske love og give python en hånd ved at give et godt estimat.
Her er den samlede kode til eksponentiel regression
Her er den samlede kode til eksponentiel regression
Klip nedenstående ud og brug i din kode husk at erstatte den røde text med dine lister
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def func(x, a, b, c):
return a * np.exp(-b * x) + c
#Her er der eksempelværdier som du skal erstatte:
x_vaerdier = np.array([0,1,2,3,4,5,6,7])
y_vaerdier = np.array([251,209,162,120,112,81,66,49])
#Erstat ovenstående tal med dine egne har du en liste slet array og skriv x_vaerdier = minliste
#print(len(x_vaerdier))
#print(len(y_vaerdier)) # tjek om de er lige lange!!!
popt, pcov = curve_fit(func, x_vaerdier, y_vaerdier, p0=[250, 0.13, 0])
a=popt[0].round(3)# popt indeholder a b og c i en liste med mange decimaler
b=popt[1].round(7)# round runder til 7 decimaler
c=popt[2].round(3)
minlegende=('y = {0} * exp(-{1}x) + {2}'.format(a,b,c))
halveringstid=(np.log(2)/b).round(3)
x_plot=np.linspace(0,len(x_vaerdier),100)
plt.plot(x_plot, func(x_plot,a,b,c), 'r-',label=minlegende)#plotte grafen
plt.plot(x_vaerdier,y_vaerdier,'ko',label='datapunkter') #plotte punkter
plt.title("x vs y = a*exp(-bx) +c")
plt.xlabel('x-værdier')
plt.ylabel('y-værdier')
plt.legend()
plt.show()
print('Regressionsligningen bliver: y={0}exp(-{1}x)+{2}'.format(a,b,c))
print('halveringstid = {0} s'.format(halveringstid))
Outputtet fra print ser sådan her ud
Hvordan virker scipy optimize?
"""
popt er en liste som indeholder a b c fra y=a*exp(-b*x)+c
pcov er covariansmatricen som vi ikke skal bruge.
curve_fit er funktionen som gør arbejdet.
1) Den fitter funktionen "func" som er e-funtionen defineret længere oppe
2) Som input bruger den x_vaerdier og y_vaerdier samt
3) en startguess for fittefunktionen (her defineret sim "func". Siden der er abc tre værdier som
definerer funktionen skal man give 3 startværdier.
4) kik på funktionen .. a er 250 .. b er ca. 0,13 og c ved vi ikke gætter på 0
"""
Page updated
Google Sites
Report abuse