Skip to content
Snippets Groups Projects
Commit 005534c9 authored by Francesco Santanastasio's avatar Francesco Santanastasio
Browse files

updated 2024

parent 67f6cfba
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
NonLinearFitExample.ipynb
+ 0
− 260
View file @ 005534c9
%% Cell type:code id:d5a34fcc tags:
``` python
# plots will be shown inline
%matplotlib inline
import matplotlib.pyplot as plt
import numpy
from numpy import sqrt,floor
import numpy as np
import pandas as pd
import random
import math
# libreria locale
import my_lib_santanastasio as my
##NOTA IMPORTANTE: se cambi qualcosa in my_lib_santanastasio
# devi fare Kernel --> Restart prima di rigirare il codice
# altrimenti i cambiamenti non saranno applicati.
from scipy import stats
import scipy.integrate as integrate
from scipy import optimize
# set global random seed
random_Gseed = 1112
np.random.seed(random_Gseed)
```
%% Cell type:code id:75ae4577 tags:
``` python
# change "virgola" (comma) with "punto" (point) in the file
myfilename = "termometro_hot_to_cold_1_60C.txt"
myfilenamemod = myfilename.split(".")[0]+str("_mod.txt")
#print (myfilenamemod)
reading_file = open(myfilename, "r", encoding="utf8", errors='ignore')
new_file_content = ""
for i, line in enumerate(reading_file):
#print (i, line)
stripped_line = line.strip()
new_line = stripped_line.replace(",", ".")
if i!=0: # drop first line of file
new_file_content += new_line +"\n"
reading_file.close()
writing_file = open(myfilenamemod, "w")
writing_file.write(new_file_content)
writing_file.close()
```
%% Cell type:code id:986c2e62 tags:
``` python
df = pd.read_csv(myfilenamemod,header=0,sep='\t') #il separatore in questo caso era un "tab" (\t)
df
```
%% Output
Tempo ( s ) Temperatura ( C )
0 0.0 60.81
1 0.1 60.70
2 0.2 60.70
3 0.3 60.70
4 0.4 60.70
.. ... ...
955 95.5 20.38
956 95.6 20.38
957 95.7 20.38
958 95.8 20.38
959 95.9 20.42
[960 rows x 2 columns]
%% Cell type:code id:fcb1d497 tags:
``` python
# Rinominare colonne
df = df.rename(columns={"Tempo ( s )": "time", "Temperatura ( C )": "temp"})
df
```
%% Output
time temp
0 0.0 60.81
1 0.1 60.70
2 0.2 60.70
3 0.3 60.70
4 0.4 60.70
.. ... ...
955 95.5 20.38
956 95.6 20.38
957 95.7 20.38
958 95.8 20.38
959 95.9 20.42
[960 rows x 2 columns]
%% Cell type:code id:a290b012 tags:
``` python
# Grafico di tutti i dati
df.plot(x="time",y="temp",linestyle="None",marker=".")
```
%% Output
<AxesSubplot:xlabel='time'>
%% Cell type:code id:f3822031 tags:
``` python
# Ripulisco i dati
dfmod = df.drop( df[ (df.time <10) | (df.time > 30) ].index )
dfmod.plot(x="time",y="temp",linestyle="None",marker=".")
```
%% Output
<AxesSubplot:xlabel='time'>
%% Cell type:code id:11abfefd tags:
``` python
t0 = dfmod.time.to_numpy()[0]
T0 = dfmod.temp.to_numpy()[0]
print (T0)
time = dfmod.time.to_numpy()-t0
temp = dfmod.temp.to_numpy()
utemp = np.repeat(0.2,len(temp)) #incertezze sulle temperature
plt.errorbar(time,temp,yerr=utemp,linestyle="None",marker=".")
plt.xlabel("Time [s]")
plt.ylabel("Temperature [C]")
```
%% Output
57.38
Text(0, 0.5, 'Temperature [C]')
%% Cell type:markdown id:a27f877e tags:
# Fit con il metodo dei minimi quadrati
%% Cell type:markdown id:400abe59 tags:
## Usando il metodo scipy.optimize.curve_fit
%% Cell type:markdown id:21e62f33 tags:
### Esempio con 2 parametri
%% Cell type:code id:fa78582c tags:
``` python
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html
from scipy.optimize import curve_fit
```
%% Cell type:code id:ea328aab tags:
``` python
def func_2par(X, tau, Tm):
Yexp = Tm-(Tm-T0)*np.exp(-X/tau)
return Yexp
xdata = time
ydata = temp
sigma_ydata = utemp
popt_2par, pcov_2par = curve_fit(func_2par, xdata, ydata, sigma=sigma_ydata, absolute_sigma=True)
print (popt_2par)
print (pcov_2par)
```
%% Output
[ 4.11298595 21.54698567]
[[ 0.00012511 -0.00022676]
[-0.00022676 0.0006979 ]]
%% Cell type:code id:d7076be0 tags:
``` python
print ("Risultato del fit con metodo dei minimi quadrati (scipy.optimize.curve_fit) - 2 parameters:")
my.PrintResult("tau_fit",popt_2par[0],np.sqrt(pcov_2par[0,0]),3,"[s]")
my.PrintResult("Tm_fit",popt_2par[1],np.sqrt(pcov_2par[1,1]),3,"[C]")
print ("Cov[tau,Tm] fit", pcov_2par[0,1], "[sC]")
print ("Rho[tau,Tm] fit", pcov_2par[0,1]/(np.sqrt(pcov_2par[0,0])*np.sqrt(pcov_2par[1,1])))
```
%% Output
Risultato del fit con metodo dei minimi quadrati (scipy.optimize.curve_fit) - 2 parameters:
tau_fit = (4.113 +/- 0.011 ) [s] [0.27%]
Tm_fit = (21.547 +/- 0.026 ) [C] [0.12%]
Cov[tau,Tm] fit -0.00022676453379438262 [sC]
Rho[tau,Tm] fit -0.7674252975925171
%% Cell type:code id:12f5d9eb tags:
``` python
# Grafico del fit
x_values = np.linspace(np.min(time),np.max(time),100)
y_values = func_2par(x_values,*popt_2par)
#y_values = popt_2par[1]-(popt_2par[1]-T0)*np.exp(-x_values/popt_2par[0]) # same as above
plt.errorbar(time,temp,yerr=utemp,marker='.',linestyle = 'None',label='data')
plt.plot(x_values,y_values,linestyle='-',linewidth=3,color='red',label='fit')
plt.xlabel("Tempo [s]")
plt.ylabel("Temperatura [C]")
plt.legend()
plt.grid()
plt.yscale('linear')
```
%% Output
%% Cell type:code id:5bf4e027 tags:
``` python
## Studio dei residui
temp_atteso = func_2par(time,*popt_2par)
#temp_atteso = popt_2par[1]-(popt_2par[1]-T0)*np.exp(-time/popt_2par[0]) # same as above
d = temp - temp_atteso
d_norm = d / utemp
##print (d)
##print (d_norm)
plt.errorbar(time,d_norm,utemp/utemp,marker='.',linestyle="")
plt.ylabel("Residui normalizzati $d/\sigma_T=(T-T_{atteso})/\sigma_T$")
plt.xlabel("Tempo [s]")
plt.grid()
```
%% Output
%% Cell type:markdown id:0a43b4ff tags:
### Esempio con 3 parametri
%% Cell type:code id:a77ac1b4 tags:
``` python
def func_3par(X, tau, Tm, Tzero):
Yexp = Tm-(Tm-Tzero)*np.exp(-X/tau)
return Yexp
xdata = time
ydata = temp
sigma_ydata = utemp
popt_3par, pcov_3par = curve_fit(func_3par, xdata, ydata, sigma=sigma_ydata, absolute_sigma=True)
print (popt_3par)
print (pcov_3par)
```
%% Output
[ 4.06770623 21.58968821 57.65442079]
[[ 0.00023173 -0.00032311 -0.00069099]
[-0.00032311 0.00078166 0.00064501]
[-0.00069099 0.00064501 0.00428592]]
%% Cell type:code id:70d51fad tags:
``` python
print ("Risultato del fit con metodo dei minimi quadrati (scipy.optimize.curve_fit) - 3 parameters:")
my.PrintResult("tau_fit",popt_3par[0],np.sqrt(pcov_3par[0,0]),3,"[s]")
my.PrintResult("Tm_fit",popt_3par[1],np.sqrt(pcov_3par[1,1]),3,"[C]")
my.PrintResult("T0_fit",popt_3par[2],np.sqrt(pcov_3par[2,2]),3,"[C]")
print ("Cov[tau,Tm] fit", pcov_3par[0,1], "[sC]")
print ("Cov[tau,T0] fit", pcov_3par[0,2], "[sC]")
print ("Cov[Tm,T0] fit", pcov_3par[1,2], "[C^2]")
print ("Rho[tau,Tm] fit", pcov_3par[0,1]/(np.sqrt(pcov_3par[0,0])*np.sqrt(pcov_3par[1,1])))
print ("Rho[tau,T0] fit", pcov_3par[0,2]/(np.sqrt(pcov_3par[0,0])*np.sqrt(pcov_3par[2,2])))
print ("Rho[Tm,T0] fit", pcov_3par[1,2]/(np.sqrt(pcov_3par[1,1])*np.sqrt(pcov_3par[2,2])))
```
%% Output
Risultato del fit con metodo dei minimi quadrati (scipy.optimize.curve_fit) - 3 parameters:
tau_fit = (4.068 +/- 0.015 ) [s] [0.37%]
Tm_fit = (21.59 +/- 0.028 ) [C] [0.13%]
T0_fit = (57.654 +/- 0.065 ) [C] [0.11%]
Cov[tau,Tm] fit -0.0003231109253328368 [sC]
Cov[tau,T0] fit -0.0006909904870936007 [sC]
Cov[Tm,T0] fit 0.0006450068923079409 [C^2]
Rho[tau,Tm] fit -0.7591875202707508
Rho[tau,T0] fit -0.6933545040572168
Rho[Tm,T0] fit 0.3523991930423434
%% Cell type:markdown id:faa7d7f4 tags:
## Usando il metodo scipy.optimize.minimize
%% Cell type:markdown id:dfb79f99 tags:
### Esempio con 2 parametri
%% Cell type:code id:5a18ee87 tags:
``` python
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize
from scipy.optimize import minimize
```
%% Cell type:code id:0c5947d0 tags:
``` python
# Chi quadro con n misure indipendenti
def ChiSquareOverTwo(theta,T0,Y,X,SigmaY):
tau,Tm = theta
Yexp = Tm-(Tm-T0)*np.exp(-X/tau) # modello di fit
#print (Yexp)
return np.sum( (Y-Yexp)**2 / SigmaY**2 ) / 2
```
%% Cell type:code id:022e6cc8 tags:
``` python
first_guesses = [1,1]
ranges = ( (1,10), (1, 100) )
result = minimize(ChiSquareOverTwo,x0=first_guesses,args=(T0,ydata,xdata,sigma_ydata),method='L-BFGS-B',bounds=ranges)
```
%% Cell type:code id:cd99c43e tags:
``` python
print ("Minimization completed with status: ", result.success)
print ("Values at minimum: ", result.x)
import numdifftools as nd
H = nd.Hessian(ChiSquareOverTwo)
COV = np.linalg.inv(H(result.x,T0,ydata,xdata,sigma_ydata))
print ("")
print ("Covariance matrix (inverse of Hessian): ")
print (COV)
```
%% Output
Minimization completed with status: True
Values at minimum: [ 4.11298602 21.54698552]
Covariance matrix (inverse of Hessian):
[[ 0.00012625 -0.00022884]
[-0.00022884 0.00070165]]
%% Cell type:code id:59ebe12a tags:
``` python
print ("Risultato del fit con metodo dei minimi quadrati (scipy.optimize.minimize + numdifftools.Hessian):")
my.PrintResult("tau_fit",result.x[0],np.sqrt(COV[0,0]),3,"[s]")
my.PrintResult("Tm_fit",result.x[1],np.sqrt(COV[1,1]),3,"[C]")
print ("Cov[tau,Tm] fit", COV[0,1], "[sC]")
print ("Rho[tau,Tm] fit", COV[0,1]/(np.sqrt(COV[0,0])*np.sqrt(COV[1,1])))
```
%% Output
Risultato del fit con metodo dei minimi quadrati (scipy.optimize.minimize + numdifftools.Hessian):
tau_fit = (4.113 +/- 0.011 ) [s] [0.27%]
Tm_fit = (21.547 +/- 0.026 ) [C] [0.12%]
Cov[tau,Tm] fit -0.00022883565311786444 [sC]
Rho[tau,Tm] fit -0.7688568296304404
%% Cell type:markdown id:c4123964 tags:
### Esempio con 3 parametri
%% Cell type:code id:085c36e5 tags:
``` python
# Chi quadro con n misure indipendenti
def ChiSquareOverTwo_3par(theta,Y,X,SigmaY):
tau,Tm,T0 = theta
Yexp = Tm-(Tm-T0)*np.exp(-X/tau) # modello di fit
#print (Yexp)
return np.sum( (Y-Yexp)**2 / SigmaY**2 ) / 2
```
%% Cell type:code id:d3ef8485 tags:
``` python
first_guesses_3par = [1,1,1]
ranges_3par = ( (1,10), (1, 100) , (1,100) )
result_3par = minimize(ChiSquareOverTwo_3par,x0=first_guesses_3par,args=(ydata,xdata,sigma_ydata),method='L-BFGS-B',bounds=ranges_3par)
```
%% Cell type:code id:770fbcab tags:
``` python
print ("Minimization completed with status: ", result_3par.success)
print ("Values at minimum: ", result_3par.x)
import numdifftools as nd
H_3par = nd.Hessian(ChiSquareOverTwo_3par)
COV_3par = np.linalg.inv(H_3par(result_3par.x,ydata,xdata,sigma_ydata))
print ("")
print ("Covariance matrix (inverse of Hessian): ")
print (COV_3par)
```
%% Output
Minimization completed with status: True
Values at minimum: [ 4.0677062 21.58968844 57.65441891]
Covariance matrix (inverse of Hessian):
[[ 0.00023771 -0.00033144 -0.00070881]
[-0.00033144 0.00079327 0.00066985]
[-0.00070881 0.00066985 0.00433903]]
%% Cell type:code id:d247ec97 tags:
``` python
print ("Risultato del fit con metodo dei minimi quadrati (scipy.optimize.minimize + numdifftools.Hessian) 3 par.:")
my.PrintResult("tau_fit",result_3par.x[0],np.sqrt(COV_3par[0,0]),3,"[s]")
my.PrintResult("Tm_fit",result_3par.x[1],np.sqrt(COV_3par[1,1]),3,"[C]")
my.PrintResult("T0_fit",result_3par.x[2],np.sqrt(COV_3par[2,2]),3,"[C]")
print ("Cov[tau,Tm] fit", COV_3par[0,1], "[sC]")
print ("Cov[tau,T0] fit", COV_3par[0,2], "[sC]")
print ("Cov[Tm,T0] fit", COV_3par[1,2], "[C^2]")
print ("Rho[tau,Tm] fit", COV_3par[0,1]/(np.sqrt(COV_3par[0,0])*np.sqrt(COV_3par[1,1])))
print ("Rho[tau,T0] fit", COV_3par[0,2]/(np.sqrt(COV_3par[0,0])*np.sqrt(COV_3par[2,2])))
print ("Rho[Tm,T0] fit", COV_3par[1,2]/(np.sqrt(COV_3par[1,1])*np.sqrt(COV_3par[2,2])))
```
%% Output
Risultato del fit con metodo dei minimi quadrati (scipy.optimize.minimize + numdifftools.Hessian) 3 par.:
tau_fit = (4.068 +/- 0.015 ) [s] [0.37%]
Tm_fit = (21.59 +/- 0.028 ) [C] [0.13%]
T0_fit = (57.654 +/- 0.066 ) [C] [0.11%]
Cov[tau,Tm] fit -0.00033144281265790374 [sC]
Cov[tau,T0] fit -0.0007088058697719556 [sC]
Cov[Tm,T0] fit 0.0006698470062342376 [C^2]
Rho[tau,Tm] fit -0.7632623692140763
Rho[tau,T0] fit -0.697923473725284
Rho[Tm,T0] fit 0.3610506281921693
%% Cell type:markdown id:c7cd5a88 tags:
# Fit integrando la posterior
%% Cell type:markdown id:7cffbbf8 tags:
### Esempio con 2 parametri
%% Cell type:code id:9396f526 tags:
``` python
# Assumo prior vaga (costante) --> posterior ~ likelihood (a meno della normalizzazione)
def posterior(tau,Tm,T0,Y,X,SigmaY):
Yexp = Tm-(Tm-T0)*np.exp(-X/tau) # modello di fit
#print (Yexp)
return np.prod( stats.norm.pdf(Y, Yexp, SigmaY) )
```
%% Cell type:markdown id:7f6b8da3 tags:
## Grafico della posterior
%% Cell type:code id:1de92ec0 tags:
``` python
from matplotlib import cm, ticker
tau_values = np.linspace(4.05,4.17, 100)
Tm_values = np.linspace(21.4,21.7, 100)
tau_mesh, Tm_mesh = np.meshgrid(tau_values, Tm_values)
posterior_values = np.zeros(len(tau_values)*len(Tm_values))
idx = 0
for tauval, Tmval in np.nditer([tau_mesh,Tm_mesh]):
posterior_values[idx]=posterior(tauval,Tmval,T0,temp,time,utemp)
idx = idx + 1
posterior_mesh = posterior_values.reshape((len(Tm_values), len(tau_values)))
fig = plt.figure()
#ax1 = plt.contourf(tau_mesh, Tm_mesh, posterior_mesh, cmap=cm.jet, locator=ticker.LogLocator(), levels=100)
import matplotlib
## - log scale
#normalize = matplotlib.colors.LogNorm()
#ax1 = plt.pcolormesh(tau_mesh, Tm_mesh, posterior_mesh, cmap=cm.jet, norm=normalize)
## - linear scale
ax1 = plt.contourf(tau_mesh, Tm_mesh, posterior_mesh, cmap=cm.jet, levels=100)
cbar = plt.colorbar(ax1)
cbar.set_label('f(tau,Tm)', rotation=270, labelpad=20)
#cbar.ax.locator_params(nbins=20)
#cbar.update_ticks()
plt.xlabel("tau [s]")
plt.ylabel("Tm [C]")
plt.grid()
plt.title("Posterior probability density function")
plt.show()
```
%% Output
%% Cell type:markdown id:24e4878a tags:
## Normalizzazione della posterior
%% Cell type:code id:13e75bc5 tags:
``` python
# Range delle variabili su cui fare inferenza
#range_tau = (3.9,4.3)
#range_Tm = (21,22)
range_tau = (4.05,4.17)
range_Tm = (21.4,21.7)
# Integrale di normalizzazione della posterior
normalization, integral_error_norm = integrate.nquad(posterior,
[range_tau, #tau
range_Tm #Tm
],
args=(T0,temp,time,utemp)
)
print ("Integrale di normalizzatione (ed errore) con scipy.integrate = ", normalization, integral_error_norm)
# Integrale a mano (come esempio)
nsteps_tau = 500
nsteps_Tm = 500
tau_values_loop = np.linspace(range_tau[0],range_tau[1], nsteps_tau)
Tm_values_loop = np.linspace(range_Tm[0],range_Tm[1], nsteps_Tm)
bin_tau = (range_tau[1] - range_tau[0]) / nsteps_tau
bin_Tm = (range_Tm[1] - range_Tm[0]) / nsteps_Tm
normalization_loop = 0
for tau_val in tau_values_loop:
for Tm_val in Tm_values_loop:
base = bin_tau * bin_Tm
altezza = posterior(tau_val,Tm_val,T0,temp,time,utemp)
normalization_loop += base * altezza
print ("Integrale di normalizzatione calcolato con un semplice loop = ", normalization_loop)
# Posterior normalizzata (PDF)
def posterior_norm(tau,Tm,T0,Y,X,SigmaY):
return posterior(tau,Tm,T0,Y,X,SigmaY) / normalization
# L'argomento dell'integrale del valore atteso di tau^k
def arg_E_tau_k(tau,Tm,k,T0,Y,X,SigmaY):
return np.power(tau,k) * posterior_norm(tau,Tm,T0,Y,X,SigmaY)
# L'argomento dell'integrale del valore atteso di Tm^k
def arg_E_Tm_k(tau,Tm,k,T0,Y,X,SigmaY):
return np.power(Tm,k) * posterior_norm(tau,Tm,T0,Y,X,SigmaY)
# L'argomento dell'integrale del valore atteso di tau*tm
def arg_E_tauTm(tau,Tm,T0,Y,X,SigmaY):
return tau * Tm * posterior_norm(tau,Tm,T0,Y,X,SigmaY)
```
%% Output
Integrale di normalizzatione (ed errore) con scipy.integrate = 1.5273411094443446e+22 43176788028416.0
Integrale di normalizzatione calcolato con un semplice loop = 1.5212378762774195e+22
%% Cell type:markdown id:3bac69d0 tags:
## Calcolo dei valori attesi, dev. std. e covarianze
%% Cell type:code id:f9c75773 tags:
``` python
E_tau, error_int_E_tau = integrate.nquad(arg_E_tau_k,
[range_tau, #tau
range_Tm #Tm
],
args=(1,T0,temp,time,utemp)#k=1
)
E_tau2, error_int_E_tau2 = integrate.nquad(arg_E_tau_k,
[range_tau, #a
range_Tm #b
],
args=(2,T0,temp,time,utemp)#k=2
)
Var_tau = E_tau2 - E_tau**2
Sigma_tau = math.sqrt(Var_tau)
```
%% Cell type:code id:d4824d94 tags:
``` python
E_Tm, error_int_E_Tm = integrate.nquad(arg_E_Tm_k,
[range_tau, #a
range_Tm #b
],
args=(1,T0,temp,time,utemp)#k=1
)
E_Tm2, error_int_E_Tm2 = integrate.nquad(arg_E_Tm_k,
[range_tau, #a
range_Tm #b
],
args=(2,T0,temp,time,utemp)#k=2
)
Var_Tm = E_Tm2 - E_Tm**2
Sigma_Tm = math.sqrt(Var_Tm)
```
%% Cell type:code id:8beecf53 tags:
``` python
E_tauTm, error_int_E_tauTm = integrate.nquad(arg_E_tauTm,
[range_tau, #tau
range_Tm #Tm
],
args=(T0,temp,time,utemp)
)
Cov_tauTm = E_tauTm - E_tau*E_Tm
Rho_tauTm = Cov_tauTm / (Sigma_tau*Sigma_Tm)
```
%% Cell type:code id:810d04db tags:
``` python
print ("Risultato del fit con integrazione della posterior:")
my.PrintResult("tau_fit",E_tau,Sigma_tau,3,"[s]")
my.PrintResult("Tm_fit",E_Tm,Sigma_Tm,3,"[C]")
print ("Cov[tau,Tm] fit", Cov_tauTm, "[sC]")
print ("Rho[tau,Tm] fit", Rho_tauTm)
```
%% Output
Risultato del fit con integrazione della posterior:
tau_fit = (4.113 +/- 0.011 ) [s] [0.27%]
Tm_fit = (21.547 +/- 0.026 ) [C] [0.12%]
Cov[tau,Tm] fit -0.0002288448051785963 [sC]
Rho[tau,Tm] fit -0.7688629634197679
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment