[PYTHON] Comparaison des programmes d'adaptation

introduction

Afin d'écrire un manuel de procédure d'analyse et un exemple de programme pour l'expérience d'un étudiant de premier cycle, j'ai essayé de trouver comment ajuster les données appropriées par la méthode du carré minimum en utilisant divers programmes, je vais donc en laisser une note.

Puisque le principe est le même en ce qu'il utilise la méthode des moindres carrés, c'est un problème si vous n'obtenez pas le même résultat, mais malheureusement je savais d'avant que la réponse diffère selon le programme, mais récemment la cause et le remède J'ai trouvé le texte qui le dit, alors j'ai décidé de l'organiser moi-même. C'est le texte qui l'a déclenché. Peter Young, "Everything you wanted to know about Data Analysis and Fitting but were afraid to ask" Apparemment, cela semble être un matériel de conférence d'une université. Ici, le `` scipy.optimize.curve_fit '' de Gnuplot et Python indique que s'il y a une erreur dans les données, la valeur d'erreur attachée au paramètre résultant est incorrecte et doit être corrigée.

Dans le cas de Gnuplot, il faut le modifier comme ça.

to get correct error bars on fit parameters from gnuplot when there are error bars on the points, you have to divide gnuplot’s asymptotic standard errors by the square root of the chi-squared per degree of freedom (which gnuplot calls FIT STDFIT and, fortunately, computes correctly).

Dans le cas de scipy.optimize '' de python, c'est compliqué, curve_fit doit être modifié, mais `` moinssq n'a pas besoin d'être modifié.

I recently learned that error bars on fit parameters given by the routine curve_fit of python also have to be corrected in the same way. This is shown in two of the python scripts in appendix H. Curiously, a different python fitting routine, leastsq, gives the error bars correctly.

Cela dit que je me demande pourquoi personne n'est au courant de cela.

It is curious that I found no hits on this topic when Googling the internet.

Exemple de données

J'ai préparé un tel ensemble de données. 20 ensembles de (x, y, ey). Je vais adapter cela en ligne droite. Je vais l'enregistrer sous forme de fichier texte data.txt```.

0.0e+00 -4.569987720595017344e-01 1.526828747143463172e+00
1.0e+00 -7.106162255269843353e-01 1.402885069270964458e+00
2.0e+00 1.105159634902675325e+00 1.735638554786020915e+00
3.0e+00 -1.939878950652441869e+00 1.011014634823069747e+00
4.0e+00 3.609690931525689983e+00 1.139915698020605550e+00
5.0e+00 8.535035219721383015e-01 9.338187791237286817e-01
6.0e+00 4.770810591544029755e+00 1.321364026236713451e+00
7.0e+00 3.323982457761388787e+00 1.703973901689593173e+00
8.0e+00 3.100622722027332578e+00 1.002313080286136637e+00
9.0e+00 4.527766245564444070e+00 9.876090792441625243e-01
1.0e+01 1.990062497396323682e+00 1.355607177365929505e+00
1.1e+01 5.113013340421659336e+00 9.283045349565146598e-01
1.2e+01 4.391676777018354905e+00 1.337677147217683160e+00
1.3e+01 5.388022504497612886e+00 9.392443558621643707e-01
1.4e+01 1.134921361159764075e+01 9.232583484294124565e-01
1.5e+01 6.067025020573844074e+00 1.186258237028150475e+00
1.6e+01 1.052771612360148445e+01 1.200732350014090954e+00
1.7e+01 6.221953870216905713e+00 8.454085761899273743e-01
1.8e+01 9.628358150028700990e+00 1.442970173161927772e+00
1.9e+01 9.493784288063746857e+00 8.196526623903285236e-01

Exemple de code dans divers programmes

RACINE DU CERN

C'est ROOT. https://root.cern.ch Le code était le plus court. Dans le cas de ROOT, s'il est simple, il est préparé sans définir de fonction d'ajustement.

fit.C

{
  TGraphErrors *g = new TGraphErrors("data.txt","%lg %lg %lg");
  g->Fit("pol1");
  g->Draw("ape");
}

Le résultat ressemble à ceci.

Gnuplot

Gnuplot est également bon pour s'adapter aux fonctions, le code requis est donc très court. Je pense que le problème avec gnuplot est que la sortie par défaut n'a pas l'air trop belle.

fit.gp

set fit errorvariables

f(x) = p0 + p1*x
fit f(x) "data.txt" u 1:2:3 via p0,p1

plot "data.txt" u 1:2:3 w yerr, f(x)

print "\n ====== ERROR CORRECTED ========"
print "Chi2/NDF = ",FIT_STDFIT**2 * FIT_NDF,"/",FIT_NDF
print "  p0 = ",p0," +- ",p0_err/FIT_STDFIT
print "  p1 = ",p1," +- ",p1_err/FIT_STDFIT

# ---------Ce qui suit est l'ajustement de l'apparence de la figure--------

set term qt font "Helvetica"
set xrange [-1:20]
set rmargin 5
set tics font ",16"
set key font ",16"
set key left top
set bars fullwidth 0
set style line 1 lc "#0080FF" lw 1 pt 7 ps 1
set style line 2 lc "#FF3333" lw 2 pt 0 ps 1

set label 1 at first 1,11 sprintf("Chi2/ndf = %5.2f / %2d",FIT_STDFIT**2 * FIT_NDF,FIT_NDF) font ",18"
set label 2 at first 1,10 sprintf("p0 = %6.3f +- %7.4f",p0,p0_err/FIT_STDFIT) font ",18"
set label 3 at first 1,9  sprintf("p1 = %6.4f +- %7.5f",p1,p1_err/FIT_STDFIT) font ",18"

plot "data.txt" u 1:2:3 w yerr ls 1,\
     f(x) ls 2

L '"Erreur standard asymptotique" écrite dans la sortie standard est incorrecte et doit être corrigée. Plus précisément, divisez la valeur d'erreur par une variable appelée FIT_STDFIT, comme dans le code ci-dessus. Si vous écrivez `` set fit errorvariables '' au début, vous pouvez également récupérer la valeur d'erreur avec la variable name_err. Si vous le modifiez, la même valeur que ROOT apparaîtra.

Final set of parameters            Asymptotic Standard Error
=======================            ==========================
p0              = -1.06859         +/- 0.9578       (89.64%)
p1              = 0.566268         +/- 0.07983      (14.1%)

correlation matrix of the fit parameters:
                p0     p1     
p0              1.000 
p1             -0.884  1.000 

 ====== ERROR CORRECTED ========
Chi2/NDF = 59.1533703771407/18
  p0 = -1.06858871709936 +- 0.528376469987239
  p1 = 0.566267669300731 +- 0.0440357299923021

Une chose à garder à l'esprit est que cette erreur ne doit être corrigée que s'il y a une erreur dans les points de données (lorsque vous spécifiez trois colonnes après `` utiliser``` lors de la collecte de données avec la commande fit). Qu'as tu besoin de faire. N'effectuez pas cette correction si tous correspondent au même poids (= données sans erreur).

python scipy.optimize.curve_fit

Le chargement des données était facile avec numpy.loadtxt```.

fit_curve_fit.py

#Lire les données
import numpy as np

data = np.loadtxt("data.txt")
xx = data.T[0] 
yy = data.T[1]
ey = data.T[2]

#Définir une fonction d'ajustement ff
def ff(x,a,b):
    return a + b*x

#Ajuster et afficher les résultats
from scipy.optimized import curve_fit
import math

par, cov = curve_fit(ff,xx,yy,sigma=ey)

chi2 = np.sum(((func_pol1(xx,par[0],par[1])-yy)/ey)**2)
print("chi2 = {:7.3f}".format(chi2))
print("p0 : {:10.5f} +- {:10.5f}".format(par[0],math.sqrt(cov[0,0]/chi2*18)))
print("p1 : {:10.5f} +- {:10.5f}".format(par[1],math.sqrt(cov[1,1]/chi2*18)))

#Affichage sur le graphique
import matplotlib.pyplot as plt

x_func = np.arange(0,20,0.1)
y_func = par[0] + par[1]*x_func

plt.errorbar(xx,yy,ey,fmt="o")
plt.plot(x_func,y_func)
plt.show()

Il semble que FIT_STDFIT de Gnuplot ne soit pas fourni, donc je vais calculer le Chi2 et le NDF par moi-même et calculer l'erreur de paramètre en utilisant la composante diagonale de la matrice de covariance de sortie. Si vous le calculez correctement, vous obtiendrez la valeur correcte.

chi2 =  59.153
p0 :   -1.06859 +-    0.52838
p1 :    0.56627 +-    0.04404

python scipy.optimize.leastsq

Je n'ai jamais utilisé ça, alors https://qiita.com/yamadasuzaku/items/6d42198793651b91a1bc J'étais autorisé à me référer. C'était un peu déroutant (pas Chi ^ 2) que ce que je devais préparer était Chi, pas la fonction que je voulais adapter.

fit_leastsq.py

#Lire les données
import numpy as np

data = np.loadtxt("data.txt")
xx = data.T[0] 
yy = data.T[1]
ey = data.T[2]

#Définir Chi
from scipy.optimize import leastsq
import math

def chi(prm,x,y,ey):
    return (((prm[0]+prm[1]*x)-y)/ey)

#Préparez la valeur initiale et l'ajustement
init_val = (-0.5, 0.5)

prm, cov, info, msg, ier = leastsq(chi,init_val,args=(xx,yy,ey),full_output=True)

chi2 = np.sum((((prm[0]+prm[1]*xx) - yy)/ey)**2)

print("chi2 = {:7.3f}".format(chi2))
print("p0 : {:10.5f} +- {:10.5f}".format(prm[0],math.sqrt(cov[0,0])))
print("p1 : {:10.5f} +- {:10.5f}".format(prm[1],math.sqrt(cov[1,1])))

L'affichage du graphique est le même que ci-dessus, il est donc omis. Le résultat est le suivant. Dans le cas de moindresq, aucune modification n'est requise, de sorte que la racine carrée de la composante diagonale de la matrice de covariance de sortie peut être utilisée telle quelle.

chi2 =  59.153
p0 :   -1.06859 +-    0.52838
p1 :    0.56627 +-    0.04404

à la fin

J'avais l'habitude de penser que le résultat de gnuplot pourrait ne pas correspondre au calcul manuel, mais j'utilise généralement ROOT, donc je ne l'ai pas vérifié sérieusement, mais j'ai finalement trouvé comment le gérer, donc c'était très rafraîchissant. Personnellement, j'y suis habitué, donc ROOT est facile, mais je pensais que Gnuplot ou Python's curve_fit serait plus facile à comprendre si j'enseignais à des étudiants de premier cycle. Cependant, les deux ont le problème que l'erreur qui doit être attachée au paramètre résultant doit être corrigée.

Au fait, je pensais qu'il serait préférable d'enseigner Python au lieu de C ou Gnuplot aux étudiants de premier cycle de nos jours, alors j'ai pensé qu'il serait préférable d'enseigner Python ainsi que ma propre étude. J'ai essayé. Bien sûr, Python est bon en ce sens qu'il peut tout faire, du traitement des données à l'affichage des graphiques, mais en ce qui concerne la définition des fonctions et l'affichage des graphiques, Gnuplot est plus intuitif et se spécialise dans le dessin de graphiques. J'ai aussi senti qu'il n'y en avait qu'un. À titre d'exemple simple, il affiche presque la même chose, mais si vous comparez les deux ci-dessous, je pense que Gnuplot est plus intuitif. Python semble mieux cependant.

gnuplot

set xrange [0:10]
f(x) = sin(x)
plot f(x)

python

import numpy as np
import matplotlib.pyplot as plt
x = np.arange(0,10,0.1)
y = np.sin(x)
plt.plot(x,y)
plt.show()