[Python] Fit Power Law

[mic_1]

Mi potreste gentilmente correggere il codice?

def OnFitPowerLaw(self, e): 
        obj = e.GetEventObject()
        self.isPressed_POWERLAW = obj.GetValue()
        x = np.array(self.degr, dtype=int) 
        y = np.array(self.cnt, dtype=int) 
        logA = np.log10(x) ; logB = np.log10(y)
        fitfunc = lambda p, x: p[0] + p[1] * x
        errfunc = lambda p, x, y: (y - fitfunc(p, x))
        pinit = [1.0, -1.0]
        out = optimize.leastsq(errfunc, pinit, args=(logA, logB), full_output=1)
        pfinal = out[0]
        index = pfinal[1]
        amp = 10.0**pfinal[0]
        xdata = np.linspace(1, max(x), len(x))    # min(x)
        powerlaw = lambda x, amp, index: amp * (x**index)
        ydata = powerlaw(xdata, amp, index)
        new_xdata = []; new_ydata = []
        print('max(y) ', max(y))
        print('max(logB) ', max(logB))
        print('max(ydata) ', max(ydata))        
        # controllare ciclo for
        for i in range(len(ydata)):
           # if ydata[i] >= 1.0 :    # min(y)):  
            if float(min(y)) <= ydata[i] <= float(max(y)) : 
                new_ydata.append(ydata[i]) ; new_xdata.append(xdata[i])
            else: continue
        #print(min(new_xdata), max(new_xdata), 'np.log10(min(new_xdata)) ', np.log10(min(new_xdata)))
        #print(min(new_ydata), max(new_ydata), 'np.log10(max(new_ydata))', np.log10(max(new_ydata)))
        self.a_POWERLAW = np.around(amp, 3) # Non è preciso
        self.b_POWERLAW = np.around(index, 3) #ok

        self.corr_POWERLAW = np.corrcoef(logA, logB, rowvar=0)[0,1]   # NON TORNA
        #self.corr_POWERLAW = np.corrcoef(new_xdata, new_ydata, rowvar=0)[0,1] 
        # Corr corretto = 0.733   NON TORNA
        #print('corr ', np.corrcoef(np.log10(xx), np.log10(yy), rowvar=0)[0,1] ) # NON TORNA

        self.R_squared_POWERLAW = np.around( linregress(logA, logB).rvalue**2, 3)  #ok
        fig, ax = plt.subplots()    
        fig.set_size_inches(int(sizeFig_x), int(sizeFig_y+2))
        #ax.loglog(xdata, ydata, "r-", linewidth=2.0)
        ax.loglog(new_xdata, new_ydata, "r-", linewidth=2.0)
        ax.scatter(self.degr, self.cnt, alpha=0.0) 
        ax.xaxis.label.set_size(10)
        ax.yaxis.label.set_size(10)
        ax.set_xscale('log')
        ax.set_yscale('log')
        ax.set_xticks(self.xticks)
        ax.set_yticks(self.yticks) 
        ax.set_xticklabels(self.xticks, fontsize=8)   
        ax.set_yticklabels(self.yticks, fontsize=8)
        self.buf3 = io.BytesIO()
        fig.savefig(self.buf3, format = 'png', dpi= 72) #, bbox_inches='tight')
        self.buf3.seek(0)
        plt.close(fig) 

Ho considerato i valori positivi che rientrano nell'intervallo dei dati iniziali, ma il FITPL viene leggermente più corto. Potreste dirmi dove sbaglio? Premetto che se carico un file diverso, il Fit Power Law viene correttamente. Grazie mille



PS: Vorrei allegare il file con i dati self.degr, self.cnt ma non so come procedere.

[mic_1]

Dimenticavo....Tutto dipende dalla condizione inserita nel ciclo, o almeno credo :

for i in range(len(ydata)):
    if float(min(y)) <= ydata[i]  : #<= float(max(y)) : 
          new_ydata.append(ydata[i]) ; new_xdata.append(xdata[i])
    else: continue

Con questa condizione ottengo la seguente immagine che viene spostata sfasando l'asse delle y:



Non so proprio come risolvere il problema. Consigli? Grazie

PS: Ho poi da risolvere il valore di correlazione che non torna e il valore di "self.a_POWERLAW " che non è preciso. Mi viene 7541.582 anzichè 7544.etc

Ringrazio anticipatamente

[vict85]

Penso che manchi qualcosa nel codice. Per esempio, mancano tutti gli import. Inoltre la funzione linregress non è definita (si tratta di scipy.stats.linregress?). Inoltre non so cosa siano sizeFig_x e sizeFig_y.

[mic_1]

Ciao, le inserisco subito :

from scipy.stats import linregress
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

Dovrei aver inserito tutte le librerie che interessano questa porzione di codice.
Le funzioni

sizeFig

le puoi sostituire con qualsiasi valore di (x, y) di dimensione figura.
Si può tranquillamente omettere la parte legata ala figura è fare un semplice plt delle 2 figure.
Tralasciando la parte della figura (da plt.subplot() per interderci ) credo si possa semplificare con

plt.loglog(new_xdata, new_ydata, "r-", linewidth=2.0)
plt.scatter(self.degr, self.cnt, alpha=0.6)
plt.show()

A me interessa principalmente la parte di calcolo. Grazie

Allego in dati in questo modo:
self.degr = [ 358, 124, 106, 96, 94, 92, 90, 88, 84, 80, 76, 74, 72, 70, 68, 64, 62, 60, 58, 56, 54, 52, 50, 48, 46, 44, 42, 40, 38, 36, 34, 32, 30, 28, 26, 24, 22, 20, 18, 16, 14, 12, 10, 8, 6, 4, 2 ]

self.cnt = [ 1, 1, 1, 1, 2, 1,, 1, 1, 1, 3, 1, 1, 3, 2, 6, 3, 3, 7, 2, 1, 1, 16, 4, 5, 9, 6, 10, 10, 9, 15, 19, 35, 21, 29, 30, 49, 43, 55, 52, 66, 70, 99, 110, 186, 238, 418, 835 ]