Resultado extraño al realizar un ajuste no lineal por mínimos cuadrados a una ley de potencia

Question

Tengo un conjunto de datos (dado a continuación en mi código MATLAB) y vs x y mi objetivo final es ajustarlo a una ley de potencia $y=ax^b$ para ver qué exponente $b$ Lo entiendo. Hice un ajuste de mínimos cuadrados no lineales y como NLLSF requiere condiciones iniciales especifiqué $b_0=10$ ya que ese era el valor esperado. Me di cuenta de que después de la instalación, $b=b_0=10$ era idéntico. No cambia en absoluto. Luego hice esto para varios valores $b_0$ entre cinco y quince y cada vez $b$ no cambia. Permanece idéntico a su valor inicial, lo cual es un poco molesto porque queremos saber cuál es el exponente. $a$ por otro lado parece converger a los mismos valores independientemente de lo que $a_0$ es.

Lo siento si esto se hace demasiado largo, pero aquí presento todo lo que he hecho.

En esta primera figura, el panel superior es el gráfico de los datos. El segundo panel son los datos (en azul) en una escala logarítmica para mostrar lo bien que deberían ajustarse a una línea de potencia. El conjunto de datos es definitivamente lineal, así que no estoy intentando ajustar algo sin sentido. Sólo por diversión, hice un ajuste lineal por mínimos cuadrados en el espacio logarítmico (en rojo) y obtuve una pendiente de 7,568 y un valor de $r^2=0.998$ lo que indica un ajuste excelente.

Entonces empecé a hacer ajustes de mínimos cuadrados no lineales en los datos. Fijé $a_0=1$ porque no parece importar. Me pasa lo mismo $a$ valores pase lo que pase. Y luego varié $b_0$ entre cinco y quince años y he aquí el resumen de los resultados.

El panel superior muestra lo que $a$ sale a por varios $b_0$ . El segundo panel muestra lo que $b$ converge a para varios $b_0$ . $b_0$ incrementos en pasos de 0,1 y en todos los casos, $b$ no cambia en absoluto. Se mantiene en su valor inicial. El tercer panel muestra $r^2$ para varios $b_0$ así que al menos la primera mitad de los ajustes son buenos, con un pico en torno al $b_0=8$ .

Así que mi pregunta es, ¿qué está pasando aquí? ¿Qué tiene de "especial" este conjunto de datos? Esto me resulta demasiado extraño como para habérmelo inventado yo. Sólo dispongo de esos seis puntos de datos, no más, así que no puedo repetir el experimento ni añadir más puntos de datos. Me parece que el espacio residual en $b$ es muy plana. Lo que sea $b_0$ que elijas, el algoritmo simplemente se queda ahí y no se mueve, pero luego por una variedad de $b_0$ también consigues un buen ajuste. ¿Cómo es posible? ¿Cómo pueden todos los exponentes entre cinco y diez digamos TODOS encajar tan bien?

¿Hay alguna forma de averiguar ese exponente? ¿Estoy condenado a hacer un ajuste lineal en el espacio logarítmico y tomar la pendiente como exponente? ¿Tiene algo que ver con el pequeño número de (seis) puntos de datos aquí? Al igual que intentar ajustar un polinomio de séptimo grado a seis puntos dará como resultado que muchos polinomios se ajusten bien. ¿Es común este fenómeno? ¿Está bien conocido o estudiado? Cualquier buena referencia también será apreciada. Gracias a todos de antemano.

También estoy pegando mi código en MATLAB en caso de que alguien quiere ver exactamente lo que estoy haciendo.

close all
clc
clear all

% The data hardcoded
x = 4:0.5:6.5;
y = 1.0e-04 * [0.011929047861902, 0.026053683604530, 0.057223759612162,...
    0.117413370572612, 0.242357468772138, 0.462327165034928];

%The data plot
figure('units','normalized','position',[0 0 1 1])

% Shows the data as it is
subplot(2,1,1)
plot(x,y,'-o')
title('The Data y vs. x to be fitted')
xlabel('x');ylabel('y')

% Shows the data on a log log scale to see how well it should fit a power
% law
subplot(2,1,2)
loglog(x,y,'-o')
hold on
title('The Data y vs. x to be fitted on a log-log scale')
xlabel('x');ylabel('y')

% LLSF in log-log space to see how well it should fit a power law
ft = fittype('poly1');
opts = fitoptions(ft);
[fitresult, gof] = fit( log(x'), log(y'), ft)
loglog(x,exp(fitresult.p2)*x.^fitresult.p1,'-r')
hold off

% Will hold the results for various initial b values
fitresultarray=[];
binit = 5:0.1:15;

% NLLSF for various initial b-values are performed and saved
for counter = 1:length(binit)

ft = fittype( 'power1' );
opts = fitoptions( ft );
opts.Display = 'Off';
opts.Lower = [-Inf -Inf];
opts.StartPoint = [1 binit(counter)];
opts.Upper = [Inf Inf];

[fitresult, gof] = fit(x',y',ft,opts);

fitresultarray = [fitresultarray; fitresult.a fitresult.b gof.rsquare];

end

%fitresultarray

% The results of the fits are plotted
figure('units','normalized','position',[0 0 1 1])
subplot(3,1,1)
plot(binit,fitresultarray(:,1),'-o')
title('Nonlinear Least Squares Fitting y vs. x with Fitting Form y=ax^b') 
xlabel('b initial');ylabel('a');
subplot(3,1,2)
plot(binit,fitresultarray(:,2),'-o')
xlabel('b initial');ylabel('b');
subplot(3,1,3)
plot(binit,fitresultarray(:,3),'-o')
xlabel('b initial');ylabel('Goodness of fit r^2');

Answer 1

Ahora que has explorado un par de formas de ajustar una ley de potencia, ¿por qué no pruebas con la forma recomendada . Para el caso continuo, el MLE para el parámetro de escala es:

$$\hat \alpha = 1 + n \left[ \sum_{i=1}^n ln \frac{x_i}{x_{min}} \right]$$

Uno de los autores ofrece consejos adicionales (algo ácidos) en su blog tiene una especie de manía con las leyes del poder y su mal uso:

Abusar de la regresión lineal hace llorar al bebé Gauss. Ajustar una línea al diagrama logarítmico por mínimos cuadrados es una mala idea. Por lo general, ni siquiera da una distribución de probabilidad, e incluso si los datos siguen una distribución de ley de potencia, da una mala estimación de los parámetros. No puedes usar las estimaciones de error que te da tu software de regresión, porque esas fórmulas incorporan suposiciones que contradicen directamente la idea de que estás viendo muestras de una ley de potencia. Y no, no se puede afirmar que porque la línea "explica" (en realidad, describe) una gran parte de la varianza que debe tener una ley de potencia, porque se puede obtener una muy alta R ^ 2 de otras distribuciones (que la prueba no tiene "poder"). Y esto sin entrar en los errores adicionales causados por tratar de ajustar una línea a histogramas binned.

Es cierto que Pareto se dedicó a trazar líneas en gráficos logarítmicos en su día, cuando inició todo este asunto de las leyes de potencia, pero "ese día" era la década de 1890. Hay un momento y un lugar para ser de la vieja escuela; este no lo es.

TL;DR Si se trata de una distribución, utilice la MLE. Compruebe la bondad del ajuste y compárelo con modelos alternativos.

Answer 2

Sería bastante inusual encontrar datos de ley de potencia para los que la varianza sea realmente constante a través de $x$ valores. Sin embargo, tomemos esa varianza constante como dada y procedamos con el análisis como un problema de mínimos cuadrados.

La forma más sencilla de ajustar este modelo concreto es utilizar GLM, no NLS. No se necesitan valores iniciales.

No tengo Matlab en esta máquina, pero en R:

x <- seq(4,6.5,0.5)
y <- 1.0e-04 * c(0.011929047861902, 0.026053683604530, 0.057223759612162,
    0.117413370572612, 0.242357468772138, 0.462327165034928)

 powfit <- glm(y~log(x),family=gaussian(link=log))
 summary(powfit)

[snip]

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -25.05254    0.16275 -153.93 1.07e-08 ***
log(x)        8.05095    0.08824   91.24 8.65e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 6.861052e-14)

    Null deviance: 1.5012e-09  on 5  degrees of freedom
Residual deviance: 2.2273e-13  on 4  degrees of freedom
AIC: -162.52

Number of Fisher Scoring iterations: 1

> plot(x,y,log="xy")
> lines(x,fitted(powfit))

Recomiendo encarecidamente el uso de GLM gaussianos con la función de enlace pertinente siempre que el modelo NLS se pueda moldear de esa forma (por ejemplo, para funciones exponenciales y de potencia).

No sé por qué tienes problemas con Matlab, ¿podría tener que ver con tus criterios de convergencia? No hay ningún problema con los datos, pero te sugiero que empieces con los coeficientes de regresión lineal:

> lm(log(y)~log(x))

Call:
lm(formula = log(y) ~ log(x))

Coefficients:
(Intercept)       log(x)  
    -24.203        7.568  

Por tanto, inicie sus parámetros en $\exp(-24.2)$ y $7.57$ y ver cómo funciona. O mejor aún, hazlo de la forma más sencilla y olvídate por completo de los valores iniciales.

Aquí están mis resultados en R utilizando nls primero, empezando por $b$ al valor del ajuste logarítmico lineal:

nls(y~a*x^b,start=list(a=exp(-24.2),b=7.57))
Nonlinear regression model
  model: y ~ a * x^b
   data: parent.frame()
        a         b 
1.290e-11 8.063e+00 
 residual sum-of-squares: 2.215e-13

Number of iterations to convergence: 10 
Achieved convergence tolerance: 1.237e-07

Y luego pasar a empezar en b=10 :

nls(y~a*x^b,start=list(a=exp(-24.2),b=10))
Nonlinear regression model
  model: y ~ a * x^b
   data: parent.frame()
        a         b 
1.290e-11 8.063e+00 
 residual sum-of-squares: 2.215e-13

Number of iterations to convergence: 19 
Achieved convergence tolerance: 9.267e-08

Como ves, están muy cerca.

El MLG no muy coinciden con las estimaciones de los parámetros, aunque un criterio de convergencia más estricto debería mejorarlo ( glm sólo da un paso antes de concluir que ha convergido por encima). Sin embargo, no hay mucha diferencia en el ajuste. Con un poco de manipulación, es posible encontrar un ajuste más cercano con el GLM (mejor que el nls), pero se empieza a llegar a desviaciones muy pequeñas (alrededor de $10^{-14}$ ), y me preocupa si realmente sólo estamos desplazando el error numérico hasta ese punto.

Answer 3

El problema con tu regresión no lineal es tu estimación inicial de A. Dices que fijaste el valor inicial en 1,0 "porque no parece importar". Trate de trazar una curva utilizando sus valores iniciales, y verá que está muy lejos de sus puntos de datos, tan lejos que el algoritmo de regresión no lineal está atascado y no puede averiguar si aumentar o disminuir los valores iniciales. Prueba 1e-10 como valor inicial para A y 5 para el valor inicial de B. Ahora la regresión no lineal funciona bien. Aquí están los resultados de GraphPad Prism: