Skip to content

Function implemented in R returns the wrong value

4 messages · Fernando de Souza Bastos, Duncan Murdoch, Peter Dalgaard

Original
Fernando de Souza Bastos
# 
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?

    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
Duncan Murdoch
#
On 10/12/2016 2:01 PM, Fernando de Souza Bastos wrote:
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
Fernando de Souza Bastos
# 
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.duncan at gmail.com>:

  
  


  


      

    

    

      
      

        


  

        

      Peter Dalgaard
    

    

      
      #
    

  


  

      
There are limits to how much people will do your debugging for you, but it looks like you are unpacking beta,gamma from par, but packing gamma,beta into start. Otherwise, print some of the values computed in the function and check.

-pd