R function returns wrong value

Asked

Viewed 57 times

1

the Log.Lik function below returns the value -INFwhen should return the value -5836.219. I can’t figure out the mistake, anyone has 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 function Log.Lik 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)
#Função log de verossimilhança
sum(ifelse(Y2>0,log(phi_t0)+log(Phi_t1)+log(1/sigma),log(1-Phi_mu2)))
r

1 answer

1


When replacing the lines below in the verisimilitude code, the same compiles!

Phi_mu2 <- pnorm(-mu2,log.p = TRUE)
 phi_t0 <- dnorm(term0,log = TRUE)
 phi_t1 <- dnorm(term1)
 Phi_t0 <- pnorm(term0,log.p = TRUE)
 Phi_t1 <- pnorm(term1,log.p = TRUE)
return(sum(ifelse(Y2>0,(phi_t0)+(Phi_t1)+log(1/sigma),(Phi_mu2))))

Browser other questions tagged r

You are not signed in. Login or sign up in order to post.