[R-br] Modelos não lineares....
Walmes Zeviani
walmeszeviani em gmail.com
Terça Agosto 9 14:14:34 BRT 2011
Thiago,
Já foi recomendado o não envio de arquivos anexos nas mensagens. Leia o guia
de postagem. Sugestão mais simples para conjuntos pequenos é usar dput()
como eu faço abaixo.
O problema observado já foi solucionado pelo Fernando. O modelo em questão é
não identificável da forma como foi escrito. Existem dois termos x^2 que o
torna não estimável. Outro ponto é que o modelo é linear. Veja o código
abaixo.
# lendo da área de transferência (ctrol+c no arquivo texto)
da <- read.table("clipboard", header=TRUE, dec=",")
str(da)
# ao dar o dput() gero a estrutura que copio do console e colo
# aqui para que todos tenham acesso imediato (código mínimo REPRODUZÍVEL)
dput(da)
# esse é o resultado do dput() que agora joguei para objeto
da <- structure(list(Arvore = 1:100,
DAP = c(8.3, 13.6, 13.4, 20.4,
10.1, 14.4, 13.1, 15.6, 13.4, 13.6, 14.4, 4.4, 10.4,
10.4, 11.6,
12.1, 12.4, 13.4, 10.6, 14.6, 12.2, 14.1, 10.4, 7.2,
4.4, 10.1,
14.6, 8.1, 15.1, 14.1, 11.6, 16.3, 11.4, 13.6, 11.4,
13.1, 13.4,
13.1, 11.4, 10.1, 12.4, 17.4, 10.3, 10.4, 7.4, 14.6,
15.3, 13.6,
10.4, 12.1, 10.6, 12.2, 9.4, 10.6, 15.3, 7.1, 16.8,
16.1, 19.6,
8.4, 12.4, 15.4, 12.6, 12.6, 12.6, 9.4, 10.5, 10.1,
13.4, 8.1,
11.6, 14.4, 11.8, 13.6, 15.6, 10.3, 8.6, 9.4, 9.1,
12.4, 11.4,
5.2, 11.6, 12.5, 7.5, 8.4, 9.3, 5.1, 18.4, 16.1, 6,
9.4, 12.5,
9.5, 11.2, 11.5, 6.5, 11.6, 9.5, 11.5),
HT = c(9, 8.65, 10.75, 11.87, 9.14, 9.95, 9.7, 10.5, 9,
8.55, 9.9,
5.15, 9.25, 8.78, 8.3, 8.8, 9.65, 10.93, 9.25, 9.48,
10.25, 10,
9.26, 7.47, 5.83, 8.4, 10.45, 7.67, 11.1, 10.95,
11.1, 11.7, 10.36, 11.15, 8.8,
0.05, 9.6, 11.5, 10.4, 10.65, 10.77, 10.62, 11.68,
8.45, 5.64,
11.37, 11.1, 12.47, 9.65, 9.85, 10.3, 9.86, 5.9,
9.92, 9.36,
8.24, 11.68, 12.13, 12.21, 8.2, 10, 11.16, 10.35,
8.1, 10.3,
8.87, 8.75, 8.8, 9.45, 8.75, 10.45, 11.82, 10.55,
10.15, 11.1,
10.23, 8.25, 10.3, 9.7, 9.72, 10.27, 5.9, 9.6, 10.1,
8.12, 8.55,
9.17, 5.47, 12.68, 10.5, 6.88, 9.72, 11.5, 10.18,
10.5, 10.26,
7.54, 8.7, 9.7, 8.43)),
.Names = c("Arvore", "DAP", "HT"),
class = "data.frame", row.names = c(NA, -100L))
dput(da) # isso é o dput(), use esse mecanismo para fornecer CMR à lista
# gráfico e chutes
a <- 6; b <- 0.16; c <- -0.125
plot(HT~DAP, da)
curve(1.30+(x^2/a+b*x+c*x^2), add=TRUE)
# esse modelo tem dois termos multiplicando x^2
# temos c*x^2 e temos (1/a)*x^2, portando, não é possível estimá-los
# um termo é função linear do outro
# gradiente do modelo não linear avaliados no ponto a, b, c
grad <- deriv3(~1.30+(x^2/a+b*x+c*x^2),
c("a","b","c"),
function(a, b, c, x) NULL)
grad.eval <- grad(a=a, b=b, c=c, x=da$DAP)
str(grad.eval)
F <- attr(grad.eval, "gradient")
eigen(t(F)%*%F) # último autovalor é praticamente zero
solve(t(F)%*%F) # sistema é singular
plot(F[,c("a","c")]) # linha perfeita diz que a é função de c (ou contrário)
# porém, a+c talvez seja uma função estimável
ac <- -.01; b <- 0.8
plot(HT~DAP, da)
curve(1.30+(ac*x^2+b*x), add=TRUE)
n0 <- nls(HT~(ac*DAP^2+b*DAP)+1.3, data=da,
start=list(ac=ac, b=b))
summary(n0)
# veja que esse modelo é linear e portanto pode-se usar lm()
# o problema é definir b0=1.3 na lm() que não sei se é possível
# perceba que foi necessário apenas uma interação para convergir
#-----------------------------------------------------------------------------
Tenho códigos de ajuste de modelos não lineares à dados de DAP nos meus
materiais (http://www.leg.ufpr.br/~walmes/cursoR/cursoR4.pdf)
À 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 em ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================
2011/8/9 Fernando Henrique Toledo <fernandohtoledo em gmail.com>
> Tiago
>
> Meus conhecimentos de MNL são torpes, mas as poucas vezes que os enfrentei
> esse erro apareceu algumas...
>
> Essa mensagem é um indicativo de que seu modelo pode estar
> superparametrizado, aí a superfície dos mínimos quadrados é plana e a matriz
> tem gradiente singular!
>
> Nesses casos vale a pena "esmiuçar" o modelo...
>
> Daqui a pouco o Walmes manda uma solução definitiva e elegante sobre isso!
>
> abraço,
> FH
>
> 2011/8/8 Thiago De paula protásio <depaulaprotasio em yahoo.com.br>
>
>> Pessoal;
>>
>> Novamente recorro a ajuda de vcs...
>> Estou tendo dificuldades no ajuste de um modelo não linear aos meus dados.
>>
>> Vejam mensagem de erro que está aparecendo.
>>
>> > a=6
>> > b=0.16
>> > c=-0.125
>>
>> > plot(HT~DAP,data=dados)
>> > curve(1.30+(x^2/a+b*x+c*x^2),add=T)
>>
>> > n <- nls(HT~(DAP^2/a+b*DAP+c*DAP^2)+1.3,data=dados,
>> + start=list(a=a, b=b, c=c))
>> Erro em nlsModel(formula, mf, start, wts) :
>> matriz gradiente singular nas estimativas iniciais dos parâmetros
>>
>> Em anexo segue a base de dados!
>>
>> Muito obrigado pela ajuda!
>>
>>
>> Thiago de Paula Protásio
>> Acadêmico de Engenharia Florestal
>> Universidade Federal de Lavras
>> Laboratório de Energia da Biomassa Florestal
>> Laboratório de Biomateriais
>> (035) 9183-2246
>>
>> _______________________________________________
>> R-br mailing list
>> R-br em 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.
>>
>
>
> _______________________________________________
> R-br mailing list
> R-br em 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.
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20110809/16dff436/attachment.html>
Mais detalhes sobre a lista de discussão R-br