[R] Function implemented in R returns the wrong value

[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

	[[alternative HTML version deleted]]