scipy.optimize.curve_fit avec une fonction log et echelle log

Version imprimable

1 pièce(s) jointe(s)

Bonjour,

J'essaie de faire un ajustement de ma courbe. Mes données brutes sont dans un fichier xlsx. Je les extrait à l'aide de pandas. Je veux faire deux ajustements différents parce qu'il y a un changement de comportement de Ra = 1e6. Nous savons que Ra est proportionnel à Nu ** a.
a = 0,25 pour Ra <1e6 et sinon a = 0,33.
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 import pandas as pd import numpy as np import matplotlib.pyplot as plt from math import log10 from scipy.optimize import curve_fit import lmfit data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float}) print(data) plt.xscale('log') plt.yscale('log') plt.scatter(data['Ra'].values, data['Nu_top'].values, label='Nu_top') plt.scatter(data['Ra'].values, data['Nu_bottom'].values, label='Nu_bottom') plt.errorbar(data['Ra'].values, data['Nu_top'].values , yerr=data['Ecart type top'].values, linestyle="None") plt.errorbar(data['Ra'].values, data['Nu_bottom'].values , yerr=data['Ecart type bot'].values, linestyle="None") def func(x,a): return 10**(np.log10(x)/a) """maxX = max(data['Ra'].values) minX = min(data['Ra'].values) maxY = max(data['Nu_top'].values) minY = min(data['Nu_top'].values) maxXY = max(maxX, maxY) parameterBounds = [-maxXY, maxXY]""" from lmfit import Model mod = Model(func) params = mod.make_params(a=0.25) ret = mod.fit(data['Nu_top'].head(10).values, params, x=data['Ra'].head(10).values) print(ret.fit_report()) popt, pcov = curve_fit(func, data['Ra'].head(10).values, data['Nu_top'].head(10).values, sigma=data['Ecart type top'].head(10).values, absolute_sigma=True, p0=[0.25]) plt.plot(data['Ra'].head(10).values, func(data['Ra'].head(10).values, *popt), 'r-', label='fit: a=%5.3f' % tuple(popt)) popt, pcov = curve_fit(func, data['Ra'].tail(4).values, data['Nu_top'].tail(4).values, sigma=data['Ecart type top'].tail(4).values, absolute_sigma=True, p0=[0.33]) plt.plot(data['Ra'].tail(4).values, func(data['Ra'].tail(4).values, *popt), 'b-', label='fit: a=%5.3f' % tuple(popt)) print(pcov) plt.grid plt.title("Nusselt en fonction de Ra") plt.xlabel('Ra') plt.ylabel('Nu') plt.legend() plt.show()
J'utilise donc le la fonction log10 dans python. Physiquement j'ai: logRa = a * logNu.
Ra = axe x
Nu = axe y
C'est pourquoi j'ai défini ma fonction func de cette manière.
Pièce jointe 480726
mes deux ajustements ne sont pas corrects comme vous pouvez le voir. J'ai une covariance égale à [0.00010971]. Il y a donc une erreur mais je comprend pas pourquoi ?
Voici le fichier de données: https://github.com/Suntoryy/Help/blob/master/data.xlsx

1 pièce(s) jointe(s)

log10(Ra) et log10(Nu)

J'ai changer le plot(x,y) en gros j'ai mis x=np.log10(Ra) et y=np.log10(Nu) afin de voir les coefficients 0.25 et 0.33 qui correspond a la pente des mes deux fits en temps normales.
Mais j'ai toujours un problème.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10, log
from scipy.optimize import curve_fit
import lmfit
 
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
plt.xscale('log')
plt.yscale('log')
plt.scatter(np.log10(data['Ra'].values), np.log10(data['Nu_top'].values), label='Nu_top')
plt.scatter(np.log10(data['Ra'].values), np.log10(data['Nu_bottom'].values), label='Nu_bottom')
 
plt.errorbar(np.log10(data['Ra'].values), np.log10(data['Nu_top'].values) , yerr=data['Ecart type top'].values, linestyle="None") 
plt.errorbar(np.log10(data['Ra'].values), np.log10(data['Nu_bottom'].values) , yerr=data['Ecart type bot'].values, linestyle="None")
 
def func(x,a):
	return a*x
 
maxX = max(data['Ra'].values)
minX = min(data['Ra'].values)
maxY = max(data['Nu_top'].values)
minY = min(data['Nu_top'].values)
maxXY = max(maxX, maxY)
parameterBounds = [-maxXY, maxXY]
 
from lmfit import Model
mod = Model(func)
params = mod.make_params(a=0.25)
ret = mod.fit(np.log10(data['Nu_top'].head(10).values), params, x=np.log10(data['Ra'].head(10).values))
print(ret.fit_report())
 
 
 
popt, pcov = curve_fit(func, np.log10(data['Ra'].head(10).values), np.log10(data['Nu_top'].head(10).values), sigma=data['Ecart type top'].head(10).values, absolute_sigma=True, p0=[0.25])
plt.plot(np.log10(data['Ra'].head(10).values), func(np.log10(data['Ra'].head(10).values), *popt), 'r-', label='fit: a=%5.3f' % tuple(popt))
 
popt, pcov = curve_fit(func, np.log10(data['Ra'].tail(4).values), np.log10(data['Nu_top'].tail(4).values), sigma=data['Ecart type top'].tail(4).values, absolute_sigma=True, p0=[0.33])
plt.plot(np.log10(data['Ra'].tail(4).values), func(np.log10(data['Ra'].tail(4).values), *popt), 'b-', label='fit: a=%5.3f' % tuple(popt))
 
print(pcov)
 
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('Ra')
plt.ylabel('Nu')
plt.legend()
plt.show()

Pièce jointe 480914

Si quelqu'un à une idée afin d'avoir un bon fit pour mes deux domaines ?

05/06/2019, 13h41
lg_53

Si vous voulez une échelle log vous écrivez

Code:

plt.xscale('log')

ok. Mais ca ne change pas la valeur des points qui constitue votre courbe. Si votre courbe a un point, disons de coordonnée (x,y) (ou y=f(x)), alors si vous tracez cela en log, le point garde ses coordonnées. Il n'a pas les coordonnées (log(x),log(y)) ! Ce sont simplement les axes de représentation qui sont distordus !

1 pièce(s) jointe(s)

Oui effectivement mais voici mon équation physique: Ra proportionnel à Nu**a
du coup : log(Ra)=a*log(Nu)
J'ai mes valeurs de Ra et Nu dans mon fichier data.xlsx.
Avec mon code j'ouvrir le fichier et je calcul log(Ra) et log(Nu) puis plot(log(Ra),log(Nu)) en log scale.
Je suis sensé avoir a=0.25 pour Ra<1e6 et sinon a=0.33

Voici le nouveau code:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10
from numpy import polyfit
import numpy.polynomial.polynomial as poly
 
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
 
x=np.log10(data['Ra'].values)
y1=np.log10(data['Nu_top'].values)
y2=np.log10(data['Nu_bottom'].values)
x2=np.log10(data['Ra'].head(11).values)
y4=np.log10(data['Nu_top'].head(11).values)
x3=np.log10(data['Ra'].tail(4).values)
y5=np.log10(data['Nu_top'].tail(4).values)
 
plt.xscale('log')
plt.yscale('log')
plt.scatter(x, y1, label='Nu_top')
plt.scatter(x, y2, label='Nu_bottom')
 
plt.errorbar(x, y1 , yerr=data['Ecart type top'].values, linestyle="None") 
plt.errorbar(x, y2 , yerr=data['Ecart type bot'].values, linestyle="None")
 
 
"""a=np.ones(10, dtype=np.float)
weights = np.insert(a,0,1E10)"""
 
 
 
coefs = poly.polyfit(x2, y4, 1)
print(coefs)
ffit = poly.polyval(x2, coefs)
plt.plot(x2, ffit, label='fit: b=%5.3f, a=%5.3f' % tuple(coefs))
 
absError = ffit - x2
 
SE = np.square(absError) # squared errors
MSE = np.mean(SE) # mean squared errors
RMSE = np.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (np.var(absError) / np.var(x2))
print('RMSE:', RMSE)
print('R-squared:', Rsquared)
print()
print('Predicted value at x=0:', ffit[0])
print()
 
 
coefs = poly.polyfit(x3, y5, 1)
ffit = poly.polyval(x3, coefs)
plt.plot(x3, ffit, label='fit: b=%5.3f, a=%5.3f' % tuple(coefs))
 
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('log10(Ra)')
plt.ylabel('log10(Nu)')
plt.legend()
plt.show()

Pièce jointe 480982

j'ai comme fit : ax+b et a comme coeff pour mes deux domaines.
Ma question est comment faire pour mettre b=0 ? afin de voir si j'ai mes deux coeff a qui tendent plus vers mes valeurs théoriques ?
Et est-ce normal que j'ai des barres d'erreurs aussi grande ?

05/06/2019, 15h32
lg_53
Ah, déjà c'est mieux, les droites recollent plus à la courbe !

Pour imposer b=0 (il me semble qu'on ne parle plus du problème initial du topic que vous avez ouvert là), dans le fit ax+b, je vois 2 possibilités

1) Ajouter le point 0,0 à vos données, et y mettre un gros poids dessus pour le fit. Ceci ne vous assureras pas b pile à 0, mais du moins il devrait en être très proche. Si cela est un problème la 2 eme solution fera b pile à 0

2)Utiliser curve_fit(func, xdata, ydata) ou vous avez défini func comme cela :
Code:

1 2 def func(x,a): return a*x

1 pièce(s) jointe(s)

Ok merci,

Oui avec la fonction func ça devrait marcher mais j'y arrive pas.

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from math import log10, log
from scipy.optimize import curve_fit
import lmfit
 
data=pd.read_excel('data.xlsx',sheet_name='Sheet2',index=False,dtype={'Ra': float})
print(data)
plt.xscale('log')
plt.yscale('log')
plt.scatter(np.log10(data['Ra'].values), np.log10(data['Nu_top'].values), label='Nu_top')
plt.scatter(np.log10(data['Ra'].values), np.log10(data['Nu_bottom'].values), label='Nu_bottom')
 
plt.errorbar(np.log10(data['Ra'].values), np.log10(data['Nu_top'].values) , yerr=data['Ecart type top'].values, linestyle="None") 
plt.errorbar(np.log10(data['Ra'].values), np.log10(data['Nu_bottom'].values) , yerr=data['Ecart type bot'].values, linestyle="None")
 
def func(x,a):
	return a*x-0.22914798835785583
 
def funcc(x,b):
	return b*x-0.8860566476931632
 
 
from lmfit import Model
mod = Model(func)
params = mod.make_params(a=0.25)
ret = mod.fit(np.log10(data['Nu_top'].head(10).values), params, x=np.log10(data['Ra'].head(10).values))
print(ret.fit_report())
 
 
popt, pcov = curve_fit(func, np.log10(data['Ra'].head(10).values), np.log10(data['Nu_top'].head(10).values), sigma=data['Ecart type top'].head(10).values, absolute_sigma=True, p0=[0.25])
plt.plot(np.log10(data['Ra'].head(10).values), func(np.log10(data['Ra'].head(10).values), *popt), 'r-', label='fit: a=%5.3f' % tuple(popt))
 
popt, pcov = curve_fit(funcc, np.log10(data['Ra'].tail(4).values), np.log10(data['Nu_top'].tail(4).values), sigma=data['Ecart type top'].tail(4).values, absolute_sigma=True, p0=[0.33])
plt.plot(np.log10(data['Ra'].tail(4).values), funcc(np.log10(data['Ra'].tail(4).values), *popt), 'b-', label='fit: a=%5.3f' % tuple(popt))
 
print(pcov)
 
plt.grid
plt.title("Nusselt en fonction de Ra")
plt.xlabel('log10(Ra)')
plt.ylabel('log10(Nu)')
plt.legend()
plt.show()

Pièce jointe 481018

J'obtiens ça et je sais pas pourquoi ...

06/06/2019, 09h05
lg_53
Une droite en échelle log n'est pas une droite (comme on le voit sur le graphe de votre post précédent, où on a une courbe qui grandit de moins en moins vite).

Donc à mon avis si vous enlevez la représentation log, vous avez qqch qui colle bien. C'est le 1ere truc à vérifier pour savoir si c'est un problème simplement d'affichage ou un problème dans le calcul du fit lui même.

Essayez de minimiser les différences entre les codes que vous présentez, car là ca change du tout au tout. Factoriser un peu les choses comme vous aviez fait là :
Code:

1 2 3 4 5 6 maxX = max(data['Ra'].values) minX = min(data['Ra'].values) maxY = max(data['Nu_top'].values) minY = min(data['Nu_top'].values) maxXY = max(maxX, maxY) parameterBounds = [-maxXY, maxXY]
Prenez l'habitude également de simplifier vos problèmes, afin de réduire le plus possible l'étau. Là vous avez un fit en 2 parties. Coupez vos données là où il faut, et regarder déjà simplement le fit de la 1ere partie. Ca fera un code moins lourd et vous et vous verez plus clair (et par la meme occasion, ceux qui vous aide aussi y verront plus clair). Vous ferez 2 parties, quand cela fonctionnera bien sur une seule, et pas avant, c'est inutile. Là vous avez voulu faire en 2 parties direct et ca coince ... Il faut faire par étape, et tester à chaque fois.