¿Cómo encontrar buenos valores de inicio para la función nls?

Question

No entiendo por qué no puedo tener una función nls para estos datos. Lo he intentado con muchos valores de inicio diferentes y siempre tengo el mismo error.

Esto es lo que he estado haciendo:

expFct2 = function (x, a, b,c) 
{ 
    a*(1-exp(-x/b)) + c 
} 
vec_x <- c(77.87,87.76,68.6,66.29) 
vec_y <- c(1,1,0.8,0.6) 
dt <- data.frame(vec_x=vec_x,vec_y=vec_y) 
ggplot(data = dt,aes(x = vec_x, y = vec_y)) + geom_point() + 
    geom_smooth(data=dt, method="nls", formula=y~expFct2(x, a, b, c), 
     se=F, start=list(a=1, b=75, c=-5) 

tengo siempre este error:

Error in method(formula, data = data, weights = weight, ...) : 
    singular gradient 

Answer 1

Esto se puede escribir con dos parámetros lineales (.lin1 y .lin2) y un parámetro no lineal (b) como este:

a*(1-exp(-x/b)) + c 
= (a+c) - a * exp(-x/b) 
= .lin1 + .lin2 * exp(-x/b) 

donde .lin1 = a+c y .lin2 = -a (por lo a = - .lin2 y c = .lin1 + .lin2) Esto nos permite use "plinear" que solo requiera la especificación de un valor inicial para el parámetro no lineal individual (eliminando el problema de cómo establecer los valores de inicio para los otros parámetros) y que converge a pesar de e valor inicial de b=75 estar lejos de la de la solución:

nls(y ~ cbind(1, exp(-x/b)), start = list(b = 75), alg = "plinear") 

Este es el resultado de una carrera de la que podemos ver en el tamaño de .lin2 que el problema está mal escalado:

> x <- c(77.87,87.76,68.6,66.29) 
> y <- c(1,1,0.8,0.6) 
> nls(y ~ cbind(1, exp(-x/b)), start = list(b = 75), alg = "plinear") 
Nonlinear regression model 
    model: y ~ cbind(1, exp(-x/b)) 
    data: parent.frame() 
     b  .lin1  .lin2 
3.351e+00 1.006e+00 -1.589e+08 
residual sum-of-squares: 7.909e-05 

Number of iterations to convergence: 9 
Achieved convergence tolerance: 9.887e-07 
> R.version.string 
[1] "R version 2.14.2 Patched (2012-02-29 r58660)" 
> win.version() 
[1] "Windows Vista (build 6002) Service Pack 2" 

EDITAR: se agregó la ejecución de muestra y se comenta sobre la escala.

Answer 2

Ajuste de un modelo no lineal de tres parámetros a cuatro puntos de datos va a ser moderadamente difícil en cualquier caso, aunque en este caso los datos se comportan bien. El punto n. ° 1 es que su valor inicial para su parámetro c (-5) estaba muy lejos. Dibujar una imagen de la curva correspondiente a sus parámetros de inicio (ver a continuación) lo ayudaría a comprender esto (por lo que reconocería que la curva que obtenga va desde c como mínimo hasta c+a como máximo, y el rango de sus datos es desde 0,6 a 1 ...)

Sin embargo, incluso con una mejor suposición inicial, me encontré a mí mismo con los parámetros de control (es decir, control=nls.control(maxiter=200)), seguido de más advertencias - nls no es conocido por su robustez. Así que probé el modelo SSasympOff, que implementa una versión autoarrancada de la curva que desea encajar.

start1 <- list(a=1, b=75, c=-5) 
start2 <- list(a=0.5, b=75, c=0.5) ## a better guess 

pfun <- function(params) { 
    data.frame(vec_x=60:90, 
      vec_y=do.call(expFct2,c(list(x=60:90),params))) 
} 
library(ggplot2) 
ggplot(data = dt,aes(x = vec_x, y = vec_y)) + geom_point() + 
    geom_line(data=pfun(start1))+ 
    geom_line(data=pfun(start2),colour="red")+ 
    geom_smooth(data=dt, method="nls", formula=y~SSasympOff(x, a, b, c), 
       se=FALSE) 

Mi consejo en general es que es más fácil de averiguar lo que está pasando y solucionar los problemas si encaja en nlsfuera de geom_smooth y construir la curva que desea añadir usando predict.nls ...

De manera más general, la forma de obtener buenos parámetros de inicio es comprender la geometría de la función que está ajustando y los parámetros que controlan qué aspectos de la curva. Como mencioné anteriormente, c es el valor mínimo de la curva de saturación-exponencial desplazada, a es el rango, y b es un parámetro de escala (puede ver que cuando x=b, la curva es 1-exp(-1) o aproximadamente 2/3 del recorrido desde el mínimo al máximo). O bien, un poco de álgebra y cálculo (es decir, límites) o jugar con la función curve() son buenas maneras de recopilar esta información.

Answer 3

Me cuesta encontrar una interpretación para sus parámetros: a es una pendiente, b la velocidad de convergencia, y a + c el límite, pero c por sí mismo no parece significar mucho. Después de volver a parametrizar su función, el problema desaparece.

f <- function (x, a,b,c) a + c * exp(-x/abs(b)) 
nls(y~f(x, a, b, c), data=dt, start=list(a=1, b=75, c=-5), trace=TRUE) 

Sin embargo, el valor de c se ve muy, muy alto: que es probablemente la razón por el modelo inicialmente no converger.

Nonlinear regression model 
    model: y ~ f(x, a, b, c) 
    data: dt 
     a   b   c 
1.006e+00 3.351e+00 -1.589e+08 
residual sum-of-squares: 7.909e-05 

Number of iterations to convergence: 9 
Achieved convergence tolerance: 2.232e-06 

Aquí hay otra parametrización más razonable de la misma función.

g <- function (x, a,b,c) a * (1-exp(-(x-c)/abs(b))) 
nls(y~g(x, a, b, c), data=dt, start=list(a=1, b=75, c=-5), trace=TRUE) 

Nonlinear regression model 
    model: y ~ g(x, a, b, c) 
    data: dt 
    a  b  c 
1.006 3.351 63.257 
residual sum-of-squares: 7.909e-05 

Number of iterations to convergence: 10 
Achieved convergence tolerance: 1.782e-06