[R] Function implemented in R returns the wrong value

[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

> On 10 Dec 2016, at 20:01 , Fernando de Souza Bastos <fsbmat at gmail.com> 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?
> 
>    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.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com