MÉTODO DE GAUSS-NEWTON
Boa tarde colegas, Estou tentando reproduzir o algoritmo de GAUSS-NEWTOM em linguagem r, entretanto não estou conseguindo automatizar os comandos. Na parte da automatização eu estou seguindo o script que o Walmes postou no site <http://ridiculas.wordpress.com/?s=newton> sobre o metodo de NEWTON-RAPHSON. A titilo de exemplo estou utilizando um ajuste para o modelo y=a*exp(b*x) conforme segue a baixo: CMR: x<-c(2,5,7,10,14,19,26,31,34,38,45,52,53,60,65) y<-c(54,50,45,37,35,25,20,16,18,13,8,11,8,4,6) a<-56.6646 #parâmetro obtido obtido pela linearização da função y=a*exp(b*x) b<-(-0.03797) #parâmetro obtido obtido pela linearização da função y=a*exp(b*x) # a função fx fornece as estimativas de (y~x) dados os coeficientes a e b. fx<-function(xi,yi,a,b,i,n){resp<-numeric(0)for(i in i:n){t<-(a*exp(b*xi[i]))resp<-c(resp,t) }resp} # a função dfa fornece as estimativas da primeira derivada (dfx/da ~ x) dado o coeficientes b. dfa<-function(xi,yi,b,i,n){resp<-numeric(0)for(i in i:n){ t<-exp(b*xi[i]) resp<-c(resp,t) } resp } # a função dfb fornece as estimativas da primeira derivada (dfx/db~x) dados os coeficientes a e b. dfb<-function(xi,yi,a,b,i,n){resp1<-numeric(0)for(i in i:n){ t1<-a*xi[i]*exp(b*xi[i]) resp1<-c(resp1,t1) } resp1 } # A matrix betas contem os chutes para a e b betas<-matrix(c(a,b),nc=1) # bcor fornece os parâmetros a e b corrigidos para a primeira interação. Estes parâmetros serão substituídos no lugar dos valores de a e b iniciais e o processo sera repetido. bcor<-(betas)+((solve((t(matrix(c((dfa(dia,diag,b,1,15)),(dfb(dia,diag,a,b,1,15))),nc=2))%*%(matrix(c((dfa(dia,diag,b,1,15)),(dfb(dia,diag,a,b,1,15))),nc=2)))))%*%((t(matrix(c((dfa(dia,diag,b,1,15)),(dfb(dia,diag,a,b,1,15))),nc=2))%*%(matrix(c(diag-(fx(dia,diag,a,b,1,15))),nc=1))))) # a função Sqerro fornece a soma de quadrados dos desvios Sqerro<-sum((diag-fx(dia,diag,a,b,1,15))^2) O problema é que eu não estou conseguindo usar os parâmetros corrigidos para fazer a nova interação: Ex: bcor[i+1]<-bcor[i]+betas[i] # este comando aproximaria mais vezes os parâmetros. dif(Sqerro[i+1]-Sqerro[i]) #Este seria o critério de parada quando a variação no Sqerro fosse insignificante (10^(-6) por exemplo) Agradeço desde já pela ajuda e peço desculpas se não fui muito claro. Att. Tiago.
Nesse script você encontra algo útil http://www.leg.ufpr.br/~walmes/cursoR/uesc/estimaci.R À disposição. Walmes. ========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br skype: walmeszeviani twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
Obrigado Walmes, vou ler o script. Gostaria de tirar mais uma duvida, eu poderia utilizar o método de Newton para tentar ajustar o meu modelo mesmo que uma das segundas derivadas parciais fossem 0 ? Ex: y = a*exp(b*x) dy/da = exp(b*x) d2y/da2 = 0 Date: Mon, 21 Jan 2013 11:06:17 -0200 From: walmeszeviani@gmail.com To: r-br@listas.c3sl.ufpr.br Subject: Re: [R-br] MÉTODO DE GAUSS-NEWTON Nesse script você encontra algo útil http://www.leg.ufpr.br/~walmes/cursoR/uesc/estimaci.R À disposição. Walmes. ==========================================================================Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br skype: walmeszevianitwitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmeslinux user number: 531218 ========================================================================== _______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forne�a c�digo m�nimo reproduz�vel.
Isso não vai ser problema não, pode haver zeros em alguns elementos da hessiana. Vários modelos se comportam como esse seu. À disposição. Walmes. ========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br skype: walmeszeviani twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
Obrigado Walmes consegui criar o loop com o script que vc me indicou. Date: Tue, 22 Jan 2013 16:49:16 -0200 From: walmeszeviani@gmail.com To: r-br@listas.c3sl.ufpr.br Subject: Re: [R-br] MÉTODO DE GAUSS-NEWTON Isso não vai ser problema não, pode haver zeros em alguns elementos da hessiana. Vários modelos se comportam como esse seu. À disposição. Walmes. ==========================================================================Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br skype: walmeszevianitwitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmeslinux user number: 531218 ========================================================================== _______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forne�a c�digo m�nimo reproduz�vel.
Walmes eu consegui criar o loop no entanto não estou conseguindo definir o critério de parada dif. CMR: dia<-c(2,5,7,10,14,19,26,31,34,38,45,52,53,60,65) diag<-c(54,50,45,37,35,25,20,16,18,13,8,11,8,4,6) d3 <- deriv3(~a*exp(b*x), c("a", "b"), function(x, a, b){ NULL }) bi<-c(55,-0.02) sqresi<-crossprod(diag-c(d3(dia,bi[1],bi[2]))) i <- 1 while(i < 10){ est <- d3(dia,bi[1],bi[2]) fx <- c(est) d <- attr(est, "gradient") sqresf <- crossprod(diag-fx) bf <- bi + (solve(t(d)%*%d)%*%(t(d)%*%(diag-fx))) bi <- c(bf) sqresi <- sqresf cat(paste(formatC(c(sqresf, bi), digits=6, format="f"), collapse="\t"), "\n") i <- i + 1 } E a proposito dei uma lida no post que você fala a respeito da anova para modelos não lineares. E fiquei com uma dúvida, seria possível testar o efeito dos betas como nos modelos lineares? Agradeço desde já pela ajuda. Att. Tiago. Date: Tue, 22 Jan 2013 16:49:16 -0200 From: walmeszeviani@gmail.com To: r-br@listas.c3sl.ufpr.br Subject: Re: [R-br] MÉTODO DE GAUSS-NEWTON Isso não vai ser problema não, pode haver zeros em alguns elementos da hessiana. Vários modelos se comportam como esse seu. À disposição. Walmes. ==========================================================================Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)Departamento de Estatística - Universidade Federal do Paraná fone: (+55) 41 3361 3573VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br skype: walmeszevianitwitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmeslinux user number: 531218 ========================================================================== _______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forne�a c�digo m�nimo reproduz�vel.
participantes (2)
-
Tiago Souza Marçal -
Walmes Zeviani