MARQUARDT in Non-linear Regression in R

1

was trying to adjust the following equation (G~G0*exp(-(b12)*t^-(b13))/exp(-(b12)*t0^-(b13)) by non-linear regression in R, from the following SAS script:

proc model data=five method=MARQUARDT maxiter=1000

My doubt arises in how to adjust this equation with the MARQUARDT method.

So far what I did in R was the following:

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)

However, I do not know how to specify the MARQUARDT algorithm, because the default one in R is the Gauss -newton algorithm in the nls function. Can you help me or do I have to use another bookstore?

Thank you very much !!

    
asked by marpra 09.07.2018 в 10:49
source

1 answer

0

You can use the minpack.lm library and compare the results with several methods present in 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

The result will be something like this:

# > 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
    
answered by 11.07.2018 в 21:08