Navigation

Inscrivez-vous gratuitement
pour pouvoir participer, suivre les réponses en temps réel, voter pour les messages, poser vos propres questions et recevoir la newsletter

Calcul scientifique Python Discussion :

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


Sujet :

Calcul scientifique Python

  1. #1
    Nouveau membre du Club
    scipy.optimize.curve_fit avec une fonction log et echelle log
    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 &#8203;&#8203;<1e6 et sinon a = 0,33.

    Code :Sélectionner tout -Visualiser dans une fenêtre à part
    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.

    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

  2. #2
    Nouveau membre du Club
    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 :Sélectionner tout -Visualiser dans une fenêtre à part
    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()




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

  3. #3
    Membre chevronné
    Si vous voulez une échelle log vous écrivez

    Code :Sélectionner tout -Visualiser dans une fenêtre à part
    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 !

  4. #4
    Nouveau membre du Club
    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 :Sélectionner tout -Visualiser dans une fenêtre à part
    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()




    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 ?

  5. #5
    Membre chevronné
    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 :Sélectionner tout -Visualiser dans une fenêtre à part
    1
    2
    def func(x,a):
         return a*x

  6. #6
    Nouveau membre du Club
    Ok merci,

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

    Code :Sélectionner tout -Visualiser dans une fenêtre à part
    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()




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

  7. #7
    Membre chevronné
    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 :Sélectionner tout -Visualiser dans une fenêtre à part
    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.

###raw>template_hook.ano_emploi###