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
    Membre à l'essai
    Homme Profil pro
    Étudiant
    Inscrit en
    novembre 2018
    Messages
    15
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Étudiant

    Informations forums :
    Inscription : novembre 2018
    Messages : 15
    Points : 15
    Points
    15
    Par défaut 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 ​​<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.
    Nom : Figure_1.png
Affichages : 113
Taille : 25,6 Ko
    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
    Membre à l'essai
    Homme Profil pro
    Étudiant
    Inscrit en
    novembre 2018
    Messages
    15
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Étudiant

    Informations forums :
    Inscription : novembre 2018
    Messages : 15
    Points : 15
    Points
    15
    Par défaut 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()
    Nom : Figure_1-1.png
Affichages : 50
Taille : 21,8 Ko

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

  3. #3
    Membre expérimenté

    Homme Profil pro
    Ingénieur calcul scientifique
    Inscrit en
    mars 2013
    Messages
    808
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Alpes Maritimes (Provence Alpes Côte d'Azur)

    Informations professionnelles :
    Activité : Ingénieur calcul scientifique

    Informations forums :
    Inscription : mars 2013
    Messages : 808
    Points : 1 604
    Points
    1 604
    Par défaut
    Si vous voulez une échelle log vous écrivez

    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
    Membre à l'essai
    Homme Profil pro
    Étudiant
    Inscrit en
    novembre 2018
    Messages
    15
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Étudiant

    Informations forums :
    Inscription : novembre 2018
    Messages : 15
    Points : 15
    Points
    15
    Par défaut
    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()
    Nom : Figure_1-2.png
Affichages : 48
Taille : 27,2 Ko

    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 expérimenté

    Homme Profil pro
    Ingénieur calcul scientifique
    Inscrit en
    mars 2013
    Messages
    808
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Alpes Maritimes (Provence Alpes Côte d'Azur)

    Informations professionnelles :
    Activité : Ingénieur calcul scientifique

    Informations forums :
    Inscription : mars 2013
    Messages : 808
    Points : 1 604
    Points
    1 604
    Par défaut
    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
    Membre à l'essai
    Homme Profil pro
    Étudiant
    Inscrit en
    novembre 2018
    Messages
    15
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Paris (Île de France)

    Informations professionnelles :
    Activité : Étudiant

    Informations forums :
    Inscription : novembre 2018
    Messages : 15
    Points : 15
    Points
    15
    Par défaut
    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()
    Nom : Figure_1-3.png
Affichages : 47
Taille : 24,1 Ko

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

  7. #7
    Membre expérimenté

    Homme Profil pro
    Ingénieur calcul scientifique
    Inscrit en
    mars 2013
    Messages
    808
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Alpes Maritimes (Provence Alpes Côte d'Azur)

    Informations professionnelles :
    Activité : Ingénieur calcul scientifique

    Informations forums :
    Inscription : mars 2013
    Messages : 808
    Points : 1 604
    Points
    1 604
    Par défaut
    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.

Discussions similaires

  1. etendre scipy.ndimage avec une fonction C
    Par erictronic dans le forum Calcul scientifique
    Réponses: 5
    Dernier message: 10/09/2010, 22h09
  2. Probleme avec une fonction d'ecriture de logs
    Par maxoudu328 dans le forum Windows
    Réponses: 3
    Dernier message: 09/05/2007, 01h48
  3. Thread avec une fonction membre d'une classe
    Par SteelBox dans le forum Windows
    Réponses: 6
    Dernier message: 01/03/2004, 02h15
  4. Retourner une valeur avec une fonction
    Par stephtbest dans le forum ASP
    Réponses: 4
    Dernier message: 31/10/2003, 17h37
  5. [VBA-E] avec une fonction value
    Par laas dans le forum Macros et VBA Excel
    Réponses: 3
    Dernier message: 28/11/2002, 14h22

Partager

Partager
  • Envoyer la discussion sur Viadeo
  • Envoyer la discussion sur Twitter
  • Envoyer la discussion sur Google
  • Envoyer la discussion sur Facebook
  • Envoyer la discussion sur Digg
  • Envoyer la discussion sur Delicious
  • Envoyer la discussion sur MySpace
  • Envoyer la discussion sur Yahoo