Scipy: curve_fit, polyfit, linregress, echelle log

Version imprimable

1 pièce(s) jointe(s)

Bonjour,

Je veux faire une fit de mes courbes avec une échelle log. Mes données brutes sont dans un fichier xlsx. Nous savons que Nu est proportionnel Ra**a et 0.25<a<0.33
fichier excel : https://github.com/Suntoryy/Help

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
61
62
63
64
65
66
import pandas as pd
import matplotlib.pyplot as plt
from math import log10
import numpy as np
import scipy
from scipy import stats
from scipy.optimize import curve_fit
import numpy.polynomial.polynomial as poly
 
data1=pd.read_excel('data.xlsx',sheet_name='Pr1-9', index=False)
data2=pd.read_excel('data.xlsx',sheet_name='Pr10-9', index=False)
data3=pd.read_excel('data.xlsx',sheet_name='Pr100-9', index=False)
data4=pd.read_excel('data.xlsx',sheet_name='Pr1000-9', index=False)
 
x=data1['Ra'].values
y1=data1['Nu_top'].values
y2=data2['Nu_top'].values
y3=data3['Nu_top'].values
y4=data4['Nu_top'].values
 
 
plt.scatter(x, y1, label='Nu_top Pr=1')
plt.scatter(x, y2, label='Nu_top Pr=10')
plt.scatter(x, y3, label='Nu_top Pr=100')
plt.scatter(x, y4, label='Nu_top Pr=1000')
 
 
plt.errorbar(x, y1 , yerr=data1['Ecart type top'].values, linestyle="None") 
plt.errorbar(x, y2 , yerr=data2['Ecart type top'].values, linestyle="None") 
plt.errorbar(x, y3 , yerr=data3['Ecart type top'].values, linestyle="None") 
plt.errorbar(x, y4 , yerr=data4['Ecart type top'].values, linestyle="None") 
 
def func(x,a):
	return x**a
 
 
popt, pcov = curve_fit(func, x, y1, sigma=data1['Ecart type top'].values)
plt.plot(x, func(x, *popt), 'r-', label='fit: a=%5.3f' % tuple(popt))
 
print(pcov)
 
coefs = poly.polyfit(x, y1, 1)
ffit = poly.polyval(x, coefs)
plt.plot(x, ffit,'b-', label='fit: b=%5.3f, a=%5.3f' % tuple(coefs))
 
 
lr = scipy.stats.linregress(x, y1)
plt.plot(x,lr[0]*x+lr[1],'g-', label='fit: a=%5.3f, b=%5.3f, R=%5.3f' % lr[0:3])
 
"""
lr = scipy.stats.linregress(x, y2)
plt.plot(x,lr[0]*x+lr[1], label='fit: a=%5.3f, b=%5.3f, R=%5.3f' % lr[0:3])
lr = scipy.stats.linregress(x, y3)
plt.plot(x,lr[0]*x+lr[1], label='fit: a=%5.3f, b=%5.3f, R=%5.3f' % lr[0:3])
lr = scipy.stats.linregress(x, y4)
plt.plot(x,lr[0]*x+lr[1], label='fit: a=%5.3f, b=%5.3f, R=%5.3f' % lr[0:3])"""
 
 
plt.xscale('log')
plt.yscale('log')
plt.grid
plt.title("Nusselt en fonction de Ra LVL 9")
plt.xlabel('Ra')
plt.ylabel('Nu')
plt.legend()
plt.show()

Voici ce que j'obtiens:
Pièce jointe 484143

Et je ne vois pas la solution pour faire un fit des mes données en echelle log. Si quelqu'un a une idée svp ?
Je vous en remercie !

17/06/2019, 17h01
Julien N
Salut,

J'ai comme l'impression que cette discussion est une redite de celle-ci https://www.developpez.net/forums/d1...g-echelle-log/.

J'ai un peu de mal à comprendre ce que vous souhaitez faire exactement. Corrigez moi si je me trompe, mais il me semble que vous chercher a dans y=x^a. Correct? Et vous essayez plusieurs méthodes curve_fit(), polyfit(), linregress(). Est-ce que c'est parce que vous n'êtes pas satisfait du résultat?

Attention à ne pas confondre fonction log, et échelle log. L'échelle log ce n'est juste qu'une représentation, et n'a donc rien à voir avec un problème de fitting.

Personnellement j'aborderais le problème comme suit:
- On sait que y=x^a. On cherche le paramètre a.
- On linéarise les données. log(y)=a*log(x).
- On utilise curve_fit(), ou autre, pour déterminer a.
- On trace ce que l'on veut. Soit y_fit=x^a
Il est plus facile pour de chercher les paramètres d'une loi linéaire avec une méthode linéaire.

J
17/06/2019, 17h14
Suntory

Bonjour,

Je suis d'accord avec vous. Je sais faire en linéarisant les données (log10(x) et log10(y)) et j'obtiens de bon résultat. Mais l'idéal aurait été de pouvoir afficher en échelle log ma courbe y_fit sur ce graphique.

Ceci est un autre problème:

En utilisant la fonction log10 mes barres d'erreurs soit mes écart types seront représenté par des barres de grande amplitude. Car log10(x) ou x est très petit <<1 on a un nombre très grand.
Je ne sais pas représenter mes données en log10 avec des barres d'erreurs correcte.

1 pièce(s) jointe(s)

Bonjour,

J'ai bien peur de ne toujours pas saisir ce qui vous gêne. Des fois je suis lent à la détente :)

Citation:

Envoyé par Suntory

Mais l'idéal aurait été de pouvoir afficher en échelle log ma courbe y_fit sur ce graphique.

Sur la figure exemple que vous avez fourni, vous êtes en échelle log... Et vos erreurs (errorbar()) tout autant. A partir du moment où vous faites plt.xscale('log') votre graph est en échelle log pour l'axe des abscisses.

Vous dites également que:

Citation:

Envoyé par Suntory

Je ne sais pas représenter mes données en log10 avec des barres d'erreurs correcte.

Est-ce que cela signifie que c'est lorsque vous utilisez une linéarisation que vous vous retrouvez avec des erreurs sur le log, et non pas sur les données bruts? Mais là encore je ne vois pas du tout ce qui vous gêne. Vous pouvez toujours déterminer vos paramètres sur x_log et y1_log, mais tracer les erreurs sur x et y1.

Navré de ne pouvoir vous aider plus.

[EDIT]: Pour mieux résoudre votre problème, regardons ce que j'en comprends:

Je charge vos données. Je travaille ici sur une seule série, donc je ne m'encombre pas d'indice pour y:
Code:

1 2 3 # Raw data data = pd.read_excel('data.xlsx', sheet_name='Pr1-9', index=False) x, y = data['Ra'].values, data['Nu_top'].values
Je cherche à fitter les données. En ce faisant je me suis rendu compte que y=x^a ne décrit pas très bien les données:
Code:

1 2 3 # Determine coefficients of y=x^a.10^b (or log_y=a*log_x+b) log_x, log_y = np.log10(x), np.log10(y) a, b = np.polyfit(log_x, log_y, 1)
Je trace sur un même graph (ax) les données d'origine, les barres d'erreurs, et le fit:
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 fig, ax = plt.subplots() # Plot fitting result along with raw data ax.scatter(x, y, label='Nu_top') ax.errorbar(x, y, yerr=data['Ecart type top'].values, linestyle='None') ax.plot(x, (x**a) *(10**b), label='fit: a=%0.3f, b=%0.3f' % (a, b)) ax.set_xscale('log') ax.set_yscale('log') ax.set_xlabel('x axis') ax.set_ylabel('y axis') plt.title('Figure') plt.legend() plt.show()
Notez que je préfère l'utilisation de fig et ax, mais que ce n'est en rien obligatoire.

J'otiens ceci:
Pièce jointe 484419

J

18/06/2019, 09h40
Suntory

Bonjour,

Merci pour ta réponse et pardon de ne pas avoir été très claire.
Je n'avais pas pensé a faire ce que vous avez fait :)
En fait j'avais fait le même chose que vous mais j'avais représenter le résultat en échelle log avec comme abscisse log(x) et ordonnée log(y). Je n'avais pas eu l'idée de ""revenir a arrière "" en faisant ceci :

Code:

ax.plot(x, (x**a) *(10**b), label='fit: a=%0.3f, b=%0.3f' % (a, b))

Du coup je me devais de calculer le log(erreur) afin que ce soit cohérent avec mes axes. D'où mon souci avec log10(x) qui me donnais des barres d'erreurs énorme alors que j'avais un tout petit écart type ce qui n'étais pas intuitif.

Merci !