Thank you Duncan, it really was that!
Fernando de Souza Bastos
Professor Assistente
Universidade Federal de Viçosa (UFV)
Campus UFV - Florestal
Doutorando em Estatística
Universidade Federal de Minas Gerais (UFMG)
Cel: (31) 99751-6586
2016-12-11 11:25 GMT-02:00 Duncan Murdoch <murdoch.dun...@gmail.com>:
> On 10/12/2016 2:01 PM, Fernando de Souza Bastos wrote:
>
>> The Log.lik function below returns the value '-INF' when it should return
>> the value -5836.219. I can not figure out the error, does anyone have any
>> suggestions?
>>
>
> I haven't read it carefully, but a likely problem is that you are using
> constructs like log(dnorm(x)) (e.g. log(phi_t0)) instead of dnorm(x, log =
> TRUE). The dnorm(x) value will underflow to zero, and taking the log will
> give you -Inf. Using the "log = TRUE" argument avoids the underflow.
>
> Duncan Murdoch
>
>
>> rm(list=ls())
>> library(ssmrob)
>> data(MEPS2001)
>> attach(MEPS2001)
>> n<-nrow(MEPS2001)
>> Log.lik <- function(par,X,W,y){
>> n <- length(y)
>> beta <- par[1:(ncol(X)+1)]
>> gamma <- par[(ncol(X)+2):(ncol(X)+ncol(W)+2)]
>> mu1 <- model.matrix(~X)%*%beta
>> mu2 <- model.matrix(~W)%*%gamma
>> sigma <- par[(ncol(X)+ncol(W)+3)]
>> rho <- par[(ncol(X)+ncol(W)+4)]
>> term0 <- (y-mu1)/sigma
>> term1 <- ((rho*(term0))+mu2)/sqrt(1-rho^2)
>> Phi_mu2 <- pnorm(mu2)
>> phi_t0 <- dnorm(term0)
>> phi_t1 <- dnorm(term1)
>> Phi_t0 <- pnorm(term0)
>> Phi_t1 <- pnorm(term1)
>> f <- (phi_t0*Phi_t1)/(sigma*Phi_mu2)
>> #Função log de verossimilhança
>>
>> return(sum(ifelse(Y2>0,log(phi_t0)+log(Phi_t1)+log(1/sigma),
>> log(1-Phi_mu2))))
>> }
>> y <- lnambx
>> Y2 <- dambexp
>> X <- cbind(age,female,educ,blhisp,totchr,ins)
>> W <- cbind(age,female,educ,blhisp,totchr,ins,income)
>>
>> gamma0=-0.6760544
>> gamma1=0.0879359
>> gamma2=0.6626647
>> gamma3=0.0619485
>> gamma4=-0.3639377
>> gamma5=0.7969515
>> gamma6=0.1701366
>> gamma7=0.0027078
>> beta0=5.0440623
>> beta1=0.2119747
>> beta2=0.3481427
>> beta3=0.0187158
>> beta4=-0.2185706
>> beta5=0.5399190
>> beta6=-0.0299875
>> sigma=1.2710176
>> rho=-0.1306012
>> beta=c(beta0,beta1,beta2,beta3,beta4,beta5,beta6)
>> gamma=c(gamma0,gamma1,gamma2,gamma3,gamma4,gamma5,gamma6,gamma7)
>>
>> start=c(gamma0,gamma1,gamma2,gamma3,gamma4,gamma5,gamma6,gam
>> ma7,beta0,beta1,beta2,beta3,beta4,beta5,beta6,sigma,rho)
>> Log.lik(start,X=X,W=W,y)
>>
>> If you run the codes below that are within the programming of the Log.lik
>> function they compile correctly!
>>
>> mu1 <- model.matrix(~X)%*%beta
>> mu2 <- model.matrix(~W)%*%gamma
>>
>> term0 <- (y-mu1)/sigma
>> term1 <- ((rho*(term0))+mu2)/sqrt(1-rho^2)
>> Phi_mu2 <- pnorm(mu2)
>> phi_t0 <- dnorm(term0)
>> phi_t1 <- dnorm(term1)
>> Phi_t0 <- pnorm(term0)
>> Phi_t1 <- pnorm(term1)
>> f <- (phi_t0*Phi_t1)/(sigma*Phi_mu2)
>> sum(ifelse(Y2>0,log(phi_t0)+log(Phi_t1)+log(1/sigma),log(1-Phi_mu2)))
>>
>>
>>
>>
>> Fernando Bastos
>>
>> [[alternative HTML version deleted]]
>>
>> ______________________________________________
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posti
>> ng-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>>
>
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
