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: curve_fit, polyfit, linregress, echelle log


Sujet :

Calcul scientifique Python

  1. #1
    Nouveau membre du Club
    Scipy: curve_fit, polyfit, linregress, echelle log
    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 :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
    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:


    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 !

  2. #2
    Membre éprouvé
    Salut,

    J'ai comme l'impression que cette discussion est une redite de celle-ci https://www.developpez.net/forums/d1978690/autres-langages/python-zope/calcul-scientifique/scipy-optimize-curve_fit-fonction-log-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

  3. #3
    Nouveau membre du Club
    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.

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


    J

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

###raw>template_hook.ano_emploi###