versucht, vernünftige Werte von scipy powerlaw fit zu bekommen

Question

Ich versuche, einige Daten aus einem Simulationscode, den ich ausgeführt habe, um eine Macht Gesetz Abhängigkeit herauszufinden. Wenn ich eine lineare Anpassung zeichne, passen die Daten nicht sehr gut.

Hier ist der Python-Skript ich die Daten passen bin mit:

#!/usr/bin/env python 
from scipy import optimize 
import numpy 

xdata=[ 0.00010851, 0.00021701, 0.00043403, 0.00086806, 0.00173611, 0.00347222] 
ydata=[ 29.56241016, 29.82245508, 25.33930469, 19.97075977, 12.61276074, 7.12695312] 

fitfunc = lambda p, x: p[0] + p[1] * x ** (p[2]) 
errfunc = lambda p, x, y: (y - fitfunc(p, x)) 

out,success = optimize.leastsq(errfunc, [1,-1,-0.5],args=(xdata, ydata),maxfev=3000) 

print "%g + %g*x^%g"%(out[0],out[1],out[2]) 

der Ausgang I erhalten: -71.205,3 + 71.174,5 * x^-9.79038e-05

Während auf der Schau dir die Passform an sieht ungefähr so ​​gut aus, wie du es von einer minimalsten Passform erwarten würdest, die Form der Ausgabe stört mich. Ich hatte gehofft, dass die Konstante nahe an der Stelle sein würde, wo man die Null erwarten würde (etwa 30). Und ich hatte erwartet, eine Leistungsabhängigkeit von einem größeren Bruchteil als 10^-5 zu finden.

Ich habe versucht, meine Daten zu reskalieren und mit den Parametern zu spielen, um optimize.leastsq ohne Glück zu spielen. Ist das, was ich erreichen will, möglich oder erlauben meine Daten das nicht? Die Berechnung ist teuer, daher ist es nicht trivial, mehr Datenpunkte zu erhalten.

Danke!

Answer 1

Es hilft, xdata zu rescale, so dass die Zahlen nicht alle so klein sind. Sie könnten in einer neuen Variablen xprime = 1000*x arbeiten. Dann passen xprime gegen y.

der kleinsten Quadrate werden Parameter finden q Einpassen

y = q[0] + q[1] * (xprime ** q[2]) 
    = q[0] + q[1] * ((1000*x) ** q[2]) 

So

p[0] = q[0] 
p[1] = q[1] * (1000**q[2]) 
p[2] = q[2] 

Dann y = p[0] + p[1] * (x ** p[2])

lassen Es hilft auch, die anfängliche Vermutung zu, etwas zu ändern näher an das gewünschte Ergebnis, wie als [max(ydata), -1, -0.5].

from scipy import optimize 
import numpy as np 

def fitfunc(p, x): 
    return p[0] + p[1] * (x ** p[2]) 
def errfunc(p, x, y): 
    return y - fitfunc(p, x) 

xdata=np.array([ 0.00010851, 0.00021701, 0.00043403, 0.00086806, 
       0.00173611, 0.00347222]) 
ydata=np.array([ 29.56241016, 29.82245508, 25.33930469, 19.97075977, 
       12.61276074, 7.12695312]) 

N = 5000 
xprime = xdata * N 

qout,success = optimize.leastsq(errfunc, [max(ydata),-1,-0.5], 
           args=(xprime, ydata),maxfev=3000) 

out = qout[:] 
out[0] = qout[0] 
out[1] = qout[1] * (N**qout[2]) 
out[2] = qout[2] 
print "%g + %g*x^%g"%(out[0],out[1],out[2]) 

ergibt

40,1253 + -282,949 * x^0,375555

Answer 2

Es ist viel besser, zuerst den Logarithmus zu nehmen, dann leastsquare auf diese lineare Gleichung passen verwenden, die Ihnen eine viel Besser passen. Es gibt ein großartiges Beispiel in der scipy cookbook, die ich unten an Ihren Code angepasst habe.

am besten passt wie diese sind: Amplitude = 0,8955, und index = -0,40943265484

Wie wir aus dem Diagramm (und Ihre Daten), wenn ihr ein Potenzgesetz passen wir den Amplitudenwert erwarten würde nicht sehen können nahe bei 30 sein. Wie in der Potenzgesetzgleichung f(x) == Amp * x ** index, also mit einem negativen Index: f(1) == Amp und f(0) == infinity.

from pylab import * 
from scipy import * 
from scipy import optimize 

xdata=[ 0.00010851, 0.00021701, 0.00043403, 0.00086806, 0.00173611, 0.00347222] 
ydata=[ 29.56241016, 29.82245508, 25.33930469, 19.97075977, 12.61276074, 7.12695312] 

logx = log10(xdata) 
logy = log10(ydata) 

# define our (line) fitting function 
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=(logx, logy), full_output=1) 

pfinal = out[0] 
covar = out[1] 

index = pfinal[1] 
amp = 10.0**pfinal[0] 

print 'amp:',amp, 'index', index 

powerlaw = lambda x, amp, index: amp * (x**index) 
########## 
# Plotting data 
########## 
clf() 
subplot(2, 1, 1) 
plot(xdata, powerlaw(xdata, amp, index))  # Fit 
plot(xdata, ydata)#, yerr=yerr, fmt='k.') # Data 
text(0.0020, 30, 'Ampli = %5.2f' % amp) 
text(0.0020, 25, 'Index = %5.2f' % index) 
xlabel('X') 
ylabel('Y') 

subplot(2, 1, 2) 
loglog(xdata, powerlaw(xdata, amp, index)) 
plot(xdata, ydata)#, yerr=yerr, fmt='k.') # Data 
xlabel('X (log scale)') 
ylabel('Y (log scale)') 

savefig('power_law_fit.png') 
show() 

Answer 3

Der normale Weg linear kleinsten Quadrate zu verwenden, um eine exponentielle Anpassung zu erhalten, ist zu tun, was fraxel suggests in his/her answer: eine gerade Linie passen einzuloggen (y_i).

Diese Methode hat jedoch numerische Nachteile, vor allem Empfindlichkeit (eine kleine Änderung in den Daten ergibt eine große Änderung in der Schätzung). Die bevorzugte Alternative ist ein nichtlinearer Ansatz der kleinsten Quadrate - er ist weniger empfindlich. Aber wenn Sie mit der linearen LS-Methode für unkritische Zwecke zufrieden sind, verwenden Sie diese einfach.