tratando de obtener valores razonables de scipy powerlaw fit

Question

Estoy tratando de ajustar algunos datos de un código de simulación que he estado ejecutando para descubrir una dependencia de la ley de potencia. Cuando trazado un ajuste lineal, los datos no encajan muy bien.

Aquí está el script en Python que estoy usando para ajustar los datos:

#!/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]) 

la salida que recibo es: -71205.3 + 71174.5 * x^-9.79038e-05

Mientras que en el trazar el ajuste parece tan bueno como cabría esperar de un ajuste leastsquares, la forma de la salida me molesta. Esperaba que la constante estuviera cerca de donde esperarías que fuera el cero (alrededor de 30). Y esperaba encontrar una dependencia de poder de una fracción mayor que 10^-5.

He intentado redimensionar mis datos y jugar con los parámetros para optimizar.leastsq sin suerte. ¿Lo que intento lograr es posible o mis datos simplemente no lo permiten? El cálculo es costoso, por lo que obtener más puntos de datos no es trivial.

Gracias!

Answer 1

Ayuda a reescalar xdata por lo que los números no son tan pequeños. Podría trabajar en una nueva variable xprime = 1000*x. Luego ajuste xprime versus y.

mínimos cuadrados encontrarán parámetros q ajuste

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

Así que

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

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

También ayuda a cambiar la estimación inicial a algo más cercano al resultado deseado, tal como [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]) 

produce

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

Answer 2

Es mucho mejor tomar primero el logaritmo, a continuación, utilizar leastsquare para adaptarse a esta ecuación lineal, lo que le dará una gran parte mejor ajuste. Hay un gran ejemplo en el scipy cookbook, que he adaptado a continuación para adaptarse a su código.

los mejores ajustes de este tipo son: Amplitud = 0,8955, y el índice = -0.40943265484

Como podemos ver en el gráfico (y sus datos), si es un ajuste de ley de potencia que no sería esperar que el valor de la amplitud Estar cerca de 30. Como en la ecuación de la ley de potencia f(x) == Amp * x ** index, con un índice negativo: f(1) == Amp y 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

La forma estándar para usar mínimos cuadrados lineales para obtener un ajuste exponencial es hacer lo fraxel suggests in his/her answer: ajustar una línea recta a log (y_i).

Sin embargo, este método tiene desventajas numéricas conocidas, particularmente sensibilidad (un pequeño cambio en los datos produce un gran cambio en la estimación). La alternativa preferida es utilizar un enfoque no lineal de mínimos cuadrados, es menos sensible. Pero si está satisfecho con el método LS lineal para fines no críticos, solo use eso.