[R] Function implemented in R returns the wrong value

[Duncan Murdoch]

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,gamma7,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 at 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.
>