- 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.
