La pregunta anterior lo dice todo. Básicamente, mi pregunta es para una función de ajuste genérica (podría ser arbitrariamente complicada) que no será lineal en los parámetros que estoy tratando de estimar, ¿cómo se eligen los valores iniciales para inicializar el ajuste? Estoy tratando de hacer mínimos cuadrados no lineales. ¿Hay alguna estrategia o método? ¿Se ha estudiado esto? ¿Alguna referencia? ¿Algo más que adivinar ad hoc? Específicamente, en este momento una de las formas de ajuste con las que estoy trabajando es una forma Gaussiana plus lineal con cinco parámetros que estoy tratando de estimar, como

donde (datos de abscisa) e y = log 10 (datos de ordenadas), lo que significa que en el espacio log-log mis datos se ven como una línea recta más una protuberancia que estoy aproximando por un gaussiano. No tengo ninguna teoría, nada que me guíe sobre cómo inicializar el ajuste no lineal, excepto tal vez gráficas y globos oculares como la pendiente de la línea y cuál es el centro / ancho de la protuberancia. Pero tengo más de cien de estos ajustes para hacerlo, en lugar de graficar y adivinar, preferiría algún enfoque que pueda automatizarse.$x = \log_{10}$$y = \log_{10}$

No puedo encontrar ninguna referencia, en la biblioteca o en línea. Lo único que se me ocurre es elegir aleatoriamente los valores iniciales. MATLAB ofrece elegir valores al azar de [0,1] distribuidos uniformemente. Entonces, con cada conjunto de datos, ejecuto el ajuste inicializado aleatoriamente miles de veces y luego elijo el que tiene el más alto . ¿Alguna otra (mejor) idea?$r^2$

Anexo # 1

Primero, aquí hay algunas representaciones visuales de los conjuntos de datos solo para mostrarles de qué tipo de datos estoy hablando. Estoy publicando los datos en su forma original sin ningún tipo de transformación y luego su representación visual en el espacio de registro-registro, ya que aclara algunas de las características de los datos mientras distorsiona otros. Estoy publicando una muestra de datos buenos y malos.

Cada uno de los seis paneles de cada figura muestra cuatro conjuntos de datos trazados juntos en rojo, verde, azul y cian, y cada conjunto de datos tiene exactamente 20 puntos de datos. Estoy tratando de ajustar cada uno de ellos con una línea recta más un gaussiano debido a los baches que se ven en los datos.

La primera figura es algunos de los buenos datos. La segunda figura es el gráfico log-log de los mismos buenos datos de la figura uno. La tercera cifra son algunos de los datos incorrectos. La cuarta figura es el gráfico log-log de la figura tres. Hay muchos más datos, estos son solo dos subconjuntos. La mayoría de los datos (aproximadamente 3/4) son buenos, similares a los buenos datos que mostré aquí.

Ahora, algunos comentarios, por favor tengan paciencia conmigo ya que esto podría ser largo, pero creo que todos estos detalles son necesarios. Trataré de ser lo más conciso posible.

Originalmente esperaba una ley de potencia simple (que significa línea recta en el espacio de registro-registro). Cuando tracé todo en el espacio log-log, vi el inesperado aumento de alrededor de 4.8 mHz. La protuberancia se investigó a fondo y se descubrió en otros trabajos, así que no es que nos equivoquemos. Está físicamente allí y otras obras publicadas mencionan esto también. Entonces, simplemente agregué un término gaussiano a mi forma lineal. Tenga en cuenta que este ajuste debía hacerse en el espacio de registro-registro (de ahí mis dos preguntas, incluida esta).

Ahora, después de leer la respuesta de Stumpy Joe Pete a otra pregunta mía (no relacionada con estos datos en absoluto) y leer esto y esto y las referencias allí (cosas de Clauset), me doy cuenta de que no debería encajar en el log-log espacio. Así que ahora quiero hacer todo en un espacio pretransformado.

Pregunta 1: Mirando los buenos datos, sigo pensando que un lineal más un gaussiano en el espacio pretransformado sigue siendo una buena forma. Me encantaría saber de otros que tienen más experiencia en datos de lo que piensan. ¿Gaussiano + lineal es razonable? ¿Debería hacer solo un gaussiano? ¿O una forma completamente diferente?

Preguntas 2: Cualquiera que sea la respuesta a la pregunta 1, aún necesitaría (muy probablemente) un ajuste de mínimos cuadrados no lineales, por lo que aún necesito ayuda con la inicialización.

Los datos donde vemos dos conjuntos, preferimos capturar el primer golpe a alrededor de 4-5 mHz. Por lo tanto, no quiero agregar más términos gaussianos y nuestro término gaussiano debería centrarse en la primera protuberancia, que casi siempre es la protuberancia más grande. Queremos "más precisión" entre 0.8mHz y alrededor de 5mHz. No nos importan demasiado las frecuencias más altas, pero tampoco queremos ignorarlas por completo. ¿Entonces quizás algún tipo de pesaje? ¿O B se puede inicializar alrededor de 4.8mHz siempre?

$f$$L$

• $f$
• $L$
• $A$$A>0$$A$
• $B$
• $C$$C$$-C$
• $D$
• $E$$L$$E$$L$$f=0$

$E$$E$$f=0$

$L$

Preguntas 3: ¿Qué piensan ustedes extrapolando de esta manera en este caso? Cualquier pros / contras? ¿Alguna otra idea para la extrapolación? Una vez más, solo nos interesan las frecuencias más bajas, por lo que extrapolamos entre 0 y 1mHz ... a veces frecuencias muy pequeñas, cercanas a cero. Sé que esta publicación ya está llena. Hice esta pregunta aquí porque las respuestas podrían estar relacionadas, pero si ustedes lo prefieren, puedo separar esta pregunta y hacer otra más tarde.

Por último, aquí hay dos conjuntos de datos de muestra a pedido.

0.813010000000000   0.091178000000000   0.012728000000000
1.626000000000000   0.103120000000000   0.019204000000000
2.439000000000000   0.114060000000000   0.063494000000000
3.252000000000000   0.123130000000000   0.071107000000000
4.065000000000000   0.128540000000000   0.073293000000000
4.878000000000000   0.137040000000000   0.074329000000000
5.691100000000000   0.124660000000000   0.071992000000000
6.504099999999999   0.104480000000000   0.071463000000000
7.317100000000000   0.088040000000000   0.070336000000000
8.130099999999999   0.080532000000000   0.036453000000000
8.943100000000001   0.070902000000000   0.024649000000000
9.756100000000000   0.061444000000000   0.024397000000000
10.569000000000001   0.056583000000000   0.025222000000000
11.382000000000000   0.052836000000000   0.024576000000000
12.194999999999999   0.048727000000000   0.026598000000000
13.008000000000001   0.045870000000000   0.029321000000000
13.821000000000000   0.041454000000000   0.067300000000000
14.633999999999999   0.039596000000000   0.081800000000000
15.447000000000001   0.038365000000000   0.076443000000000
16.260000000000002   0.036425000000000   0.075912000000000

La primera columna son las frecuencias en mHz, idénticas en cada conjunto de datos. La segunda columna es un buen conjunto de datos (datos buenos, figuras uno y dos, panel 5, marcador rojo) y la tercera columna es un conjunto de datos incorrectos (datos malos, figuras tres y cuatro, panel 5, marcador rojo).

Espero que esto sea suficiente para estimular una discusión más iluminada. Gracias a todos.

Si hubiera una estrategia que fuera a la vez buena y general , una que siempre funcionara, ya se implementaría en todos los programas de mínimos cuadrados no lineales y los valores iniciales no serían un problema.

Para muchos problemas específicos o familias de problemas, existen algunos enfoques bastante buenos para los valores iniciales; algunos paquetes vienen con buenos cálculos de valores iniciales para modelos no lineales específicos o con enfoques más generales que a menudo funcionan pero que pueden necesitar ayuda con funciones más específicas o entrada directa de valores iniciales.

Explorar el espacio es necesario en algunas situaciones, pero creo que es probable que su situación sea tal que las estrategias más específicas valen la pena, pero diseñar una buena requiere mucho conocimiento de dominio que es poco probable que poseamos.

$x$

$y$$x$

$A$

Algunos datos de muestra ayudarían, casos típicos y difíciles, si puede.

Editar: Aquí hay un ejemplo de cómo puede hacerlo bastante bien si el problema no es demasiado ruidoso:

Aquí hay algunos datos que se generan a partir de su modelo (los valores de población son A = 1.9947, B = 10, C = 2.828, D = 0.09, E = 5):

Los valores iniciales que pude estimar son
(As = 1.658, Bs = 10.001, Cs = 3.053, Ds = 0.0881, Es = 5.026)

El ajuste de ese modelo de inicio se ve así:

Los pasos fueron:

1. Ajuste una regresión de Theil para obtener una estimación aproximada de D y E
2. Resta el ajuste de la regresión de Theil
3. Use LOESS para ajustar una curva suave
4. Encuentre el pico para obtener una estimación aproximada de A, y el valor de x correspondiente al pico para obtener una estimación aproximada de B
5. Tome los ajustes LOESS cuyos valores de y son> 60% de la estimación de A como observaciones y ajuste un valor cuadrático
6. Use la cuadrática para actualizar la estimación de B y estimar C
7. De los datos originales, reste la estimación del gaussiano
8. Ajuste otra regresión de Theil a los datos ajustados para actualizar la estimación de D y E

En este caso, los valores son muy adecuados para comenzar un ajuste no lineal.

Escribí esto como Rcódigo, pero lo mismo podría hacerse en MATLAB.

Creo que cosas mejores que esto son posibles.

Si los datos son muy ruidosos, esto no funcionará del todo bien.

Edit2: este es el código que utilicé en R, si alguien está interesado:

gausslin.start <- function(x,y) {

  theilreg <- function(x,y){
    yy <- outer(y, y, "-")
    xx <- outer(x, x, "-")
    z  <- yy / xx
    slope     <- median(z[lower.tri(z)])
    intercept <- median(y - slope * x)
    cbind(intercept=intercept,slope=slope)
  }

  tr <- theilreg(x,y1)
  abline(tr,col=4)
  Ds = tr[2]
  Es = tr[1]
  yf  <- y1-Ds*x-Es
  yfl <- loess(yf~x,span=.5)

  # assumes there are enough points that the maximum there is 'close enough' to 
  #  the true maximum

  yflf   <- yfl$fitted    
  locmax <- yflf==max(yflf)
  Bs     <- x[locmax]
  As     <- yflf[locmax]

  qs     <- yflf>.6*As
  ys     <- yfl$fitted[qs]
  xs     <- x[qs]-Bs
  lf     <- lm(ys~xs+I(xs^2))
  bets   <- lf$coefficients
  Bso    <- Bs
  Bs     <-  Bso-bets[2]/bets[3]/2
  Cs     <- sqrt(-1/bets[3])
  ystart <- As*exp(-((x-Bs)/Cs)^2)+Ds*x+Es

  y1a <- y1-As*exp(-((x-Bs)/Cs)^2)
  tr  <- theilreg(x,y1a)
  Ds  <- tr[2]
  Es  <- tr[1]
  res <- data.frame(As=As, Bs=Bs, Cs=Cs, Ds=Ds, Es=Es)
  res
}

.

# population parameters: A = 1.9947 , B = 10, C = 2.828, D = 0.09, E = 5
# generate some data
set.seed(seed=3424921)
x  <- runif(50,1,30)
y  <- dnorm(x,10,2)*10+rnorm(50,0,.2)
y1 <- y+5+x*.09 # This is the data
xo <- order(x)

starts <- gausslin.start(x,y1)
ystart <- with(starts, As*exp(-((x-Bs)/Cs)^2)+Ds*x+Es)
plot(x,y1)
lines(x[xo],ystart[xo],col=2)

Existe un enfoque general para ajustar este tipo de modelos no lineales. Implica volver a parametrizar los parámetros lineales con valores de la variable dependiente, por ejemplo, el primer, último valor de frecuencia y un buen punto en el medio, digamos el sexto punto. entonces puede mantener estos parámetros fijos y resolver el parámetro no lineal en la primera fase de la minimización y luego minimizar los 5 parámetros generales.

Schnute y yo descubrimos esto alrededor de 1982 cuando ajustamos modelos de crecimiento para peces.

http://www.nrcresearchpress.com/doi/abs/10.1139/f80-172

Sin embargo, no es necesario leer este documento. Debido al hecho de que los parámetros son lineales, simplemente es necesario configurar y resolver un sistema lineal de ecuaciones de 3x3 para usar la parametrización estable del modelo.

$M$

DATA_SECTION
  init_int n
  int mid
 !! mid=6;
  init_matrix data(1,n,1,3)
  vector x(1,n)
  vector y(1,n)
 !! x=column(data,1);
 !! y=column(data,3);   //use column 3
PARAMETER_SECTION
  init_number L1(3)     //(3) means estimate in phase 3
  init_number Lmid(3)
  init_number Ln(3)

  vector L(1,3)
  init_number log_B       // estimate in phase 1
  init_number log_C(2)    // estimate in phase 2 
  matrix M(1,3,1,3);
  objective_function_value f
  sdreport_vector P(1,3)
  sdreport_number B
  sdreport_number C
  vector pred(1,n);
PROCEDURE_SECTION
  L(1)=L1;
  L(2)=Lmid;
  L(3)=Ln;
  B=exp(log_B);
  C=exp(log_C);
  M(1,1)=exp(-square((x(1)-B)/C));
  M(1,2)=x(1);
  M(1,3)=1;
  M(2,1)=exp(-square((x(mid)-B)/C));
  M(2,2)=x(mid);
  M(2,3)=1;
  M(3,1)=exp(-square((x(n)-B)/C));
  M(3,2)=x(n);
  M(3,3)=1;

  P=solve(M,L);  // solve for standard parameters 
                 // P is vector corresponding to A,D,E

  pred=P(1)*exp(-square((x-B)/C))+P(2)*x+P(3);
  if (current_phase()<4)
    f+=norm2(y-pred);
  else
    f+=0.5*n*log(norm2(y-pred))  //concentrated likelihood

$B$$C$$B$$B$$C$

Para su caso con los datos incorrectos, se ajusta con bastante facilidad y las estimaciones de parámetros (habituales) son:

         estimate    std dev
A      2.0053e-01 5.8723e-02
D      1.6537e-02 4.7684e-03
E     -1.8197e-01 7.3355e-02
B      3.0609e+00 5.0197e-01
C      5.6154e+00 9.4564e-01]

Si tiene que hacer esto muchas veces, le sugiero que use un Algoritmo Evolutivo en la función SSE como front-end para proporcionar los valores iniciales.

Por otro lado, podría usar GEOGEBRA para crear la función usando controles deslizantes para los parámetros y jugar con ellos para obtener valores iniciales.

O los valores iniciales de los datos se pueden estimar por observación.

1. D y E provienen de la pendiente e intercepción de los datos (ignorando el gaussiano)
2. A es la distancia vertical del máximo de Gauss desde la estimación de la línea Dx + E.
3. B es el valor x del máximo del gaussiano
4. C es la mitad del ancho aparente del gaussiano

Para los valores iniciales, puede hacer un ajuste de mínimos cuadrados ordinarios. Su pendiente e intersección serían los valores iniciales para D y E. El mayor residuo sería el valor inicial para A. La posición del mayor residuo sería el valor inicial para B. Quizás alguien más pueda sugerir un valor inicial para sigma.

Sin embargo, los mínimos cuadrados no lineales sin derivar ningún tipo de ecuación mecanicista del conocimiento de la materia es un negocio arriesgado, y hacer muchos ajustes por separado hace que las cosas sean aún más cuestionables. ¿Hay algún conocimiento de la materia detrás de su ecuación propuesta? ¿Hay otras variables independientes que se relacionan con las diferencias entre los 100 ajustes separados? Puede ser útil si puede incorporar esas diferencias en una sola ecuación que se ajuste a todos los datos a la vez.