Fernando de Souza Bastos
2016-Dec-10 19:01 UTC
[R] Function implemented in R returns the wrong value
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]]
Duncan Murdoch
2016-Dec-11 13:25 UTC
[R] Function implemented in R returns the wrong value
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. >
Fernando de Souza Bastos
2016-Dec-11 14:25 UTC
[R] Function implemented in R returns the wrong value
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>:> 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,gam >> ma7,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/posti >> ng-guide.html >> and provide commented, minimal, self-contained, reproducible code. >> >> >[[alternative HTML version deleted]]
peter dalgaard
2016-Dec-11 20:47 UTC
[R] Function implemented in R returns the wrong value
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