1

estaba intentando ajustar la siguiente ecuación (G~G0*exp(-(b12)*t^-(b13))/exp(-(b12)*t0^-(b13)) por regresión no lineal en R, a partir del siguiente script de SAS:

proc model data=five method=MARQUARDT maxiter=1000

Mi duda surge en como ajustar dicha ecuación con el método MARQUARDT.

Hasta ahora lo que hice en R fue lo siguiente:

library("nlstools")
nls=nls(G~G0*exp(-(b12)*t^-(b13))/exp(-(b12)*t0^-(b13)),data=ajuste,start=list(b12=3.405398, b13=0.40278))
summary(nls)
resi=nlsResiduals(nls)
plot(resi)

Sin embargo, nose como especificar el algoritmo de MARQUARDT, porque el que viene por defecto en R es el de Gauss -newton en la función nls. Me pueden ayudar o tengo que usar otra librería?

Muchas gracias!!

1

1 respuesta 1

0

Puedes utilizar la librería minpack.lm y comparar los resultados con vários métodos presentes en optimx

#preliminares, datos y librerías

library("nlstools")
library(optimx)
library(minpack.lm)

G0=10
t0=2
t=2:20
b12=3.405398
b13=0.40278
G=G0*exp(-(b12)*t^-(b13))/exp(-(b12)*t0^-(b13))+rnorm(19,0,0.2)

#Funciones

Gfun=function(x,t,G0,t0){ 
  Gc=G0*exp(-(x[1])*t^-(x[2]))/exp(-(x[1])*t0^-(x[2]))
  return(Gc)
}

Gresid <- function(parS, GObs, t,G0,t0) { # Resíduos
  GObs-Gfun(parS,t=t,G0=G0,t0=t0)
}

parStart=c(.2,.4)

ssqfun(parStart,G,t,G0,t0) # parS=parStart  GObs=G

#Levenberg-Marquardt 

nls.out <- nls.lm(par=parStart, fn = Gresid, GObs = G, t = t,G0=G0,t0=t0,
                  control = nls.lm.control(nprint=1,
                                           ftol = .Machine$double.eps,
                                           ptol = .Machine$double.eps,
                                           maxfev=10000, maxiter = 500))
summary(nls.out)

ssqfun <- function(parS, GObs, t,G0,t0) { # suma de resíduos al cuadrado
    sum((GObs-Gfun(parS,t=t,G0=G0,t0=t0))^2)
}

res<-optimx(parStart,fn=ssqfun, 
             control=list(all.methods=TRUE, save.failures=TRUE, trace=0), t=t,G0=G0,t0=t0,G=G)
res

El resultado será algo así:

# > summary(nls.out)
# 
# Parameters:
#   Estimate Std. Error t value Pr(>|t|)    
# [1,] 3.353064   0.018053  185.74   <2e-16 ***
# [2,] 0.419175   0.006355   65.96   <2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.1538 on 17 degrees of freedom
# 
# No todos optimizan -> convcode  diferente de 0
# > res
#                p1        p2         value       fevals gevals niter convcode  kkt1  kkt2 xtimes
# BFGS         3.352734 0.4192914  4.023068e-01    197     47    NA        0 FALSE  TRUE   0.00
# CG          15.433766 3.5751572  1.169458e+03    503    101    NA        1 FALSE FALSE   0.00
# Nelder-Mead  3.352300 0.4194201  4.023563e-01     97     NA    NA        0 FALSE  TRUE   0.00
# L-BFGS-B     3.352734 0.4192921  4.023069e-01     45     45    NA        0 FALSE  TRUE   0.00
# nlm          5.101472 0.1906924  3.428862e+01     NA     NA   100        1 FALSE  TRUE   0.01
# nlminb       3.353064 0.4191751  4.022989e-01     35     51    22        0  TRUE  TRUE   0.00
# spg          3.353068 0.4191737  4.022989e-01    108     NA    93        0 FALSE  TRUE   0.16
# ucminf       3.353006 0.4191950  4.022991e-01     26     26    NA        0 FALSE  TRUE   0.00
# Rcgmin             NA        NA 8.988466e+307     NA     NA    NA     9999    NA    NA   0.00
# Rvmmin             NA        NA 8.988466e+307     NA     NA    NA     9999    NA    NA   0.00
# newuoa       3.353064 0.4191753  4.022989e-01    147     NA    NA        0  TRUE  TRUE   0.00
# bobyqa       3.353064 0.4191753  4.022989e-01    245     NA    NA        0  TRUE  TRUE   0.00
# nmkb         3.353022 0.4191850  4.022997e-01    154     NA    NA        0 FALSE  TRUE   0.00
# hjkb         0.200000 0.4000000  1.200466e+04      1     NA     0     9999    NA    NA   0.00
2
  • Muchas gracias, pero no entiendo lo que significa ftol y ptol en la funcion que se llama nls.out. Me lo podrían explicar?gracias
    – marpra
    Commented el 12 jul. 2018 a las 8:32
  • Son condiciones para finalizar a optimización, definidas en help(nls.lm.control)
    – Robert
    Commented el 13 jul. 2018 a las 14:19

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.