[R-br] Diferença de resultados entre 64 e 32 bits
Rodrigo Coster
rcoster em gmail.com
Segunda Julho 16 21:04:47 BRT 2012
Caros,
programei uma rotina para estimar por máxima verossimilhança
os parâmetros de uma cópula e para ver se estava certo comparei os
resultados com o comando do pacote copula. Encontrei diferenças apenas na
5a casa decimal em diante, que considerei como sendo por causa do método
numérico utilizado. Só que, ao fazer a mesma comparação num computador 64
bits os resultados são bastante divergentes (o meu código muda o valor
estimado, enquanto o pacote copula mantem), mudando na 2a casa decimal (no
caso o parâmetro é a correlação, então a 2a casa decimal é bem importante).
Dai me bateu a seguinte dúvida: em qual confiar?
Código:
require(copula)
require(MASS)
set.seed(31415)
normCop <- function(param,data) {
n <- nrow(data)
if (length(param) != n) { param <- rep(param,nrow(data)) }
cop <- mapply(normalCopula,param=param,MoreArgs=list(dim=2))
datalist <- apply(data,1,list)
for (i in 1:n) { datalist[[i]] <- datalist[[i]][[1]] }
out <- -sum(log(mapply(dcopula,copula=cop,u=datalist)))
if (out == Inf) { out = exp(100) }
return(out)
}
a <- matrix(0,20,2)
n <- 100
Sigma <- matrix(c(10,3,3,2),2,2)
for (j in 1:20) {
data <- mvrnorm(n=n, rep(0, 2), Sigma)
data <- apply(data,2,rank)/(n+1)
fitNormCop <- function(data) {
optim(cor(data)[2],normCop,data=data, lower = 0, upper =
.9999,method="L-BFGS-B")
}
a[j,1] <- fitNormCop(data)$par # Meu
a[j,2] <- fitCopula(normalCopula(.2,2), data, method="ml")@estimate #
Pacote copula
}
E aqui um gráfico de dispersão comparando todos:
http://img411.imageshack.us/img411/7527/92588280.png
[]'s
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20120716/610a3eb5/attachment.html>
Mais detalhes sobre a lista de discussão R-br