.packageName <- "mombf"
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/dimom.r"
dimom <- function(x,V1=1,g=1,n=1,nu=1,theta0,logscale=FALSE) {

if (is.vector(x)) {
  if (missing(theta0)) theta0 <- 0
  if (missing(V1)) V1 <- 1
  qtheta <- (x-theta0)^2/(n*g*V1)
  k <- -0.5*log(n*g*V1) - lgamma(nu/2)
  p1 <- 1
} else {
  if (missing(theta0)) theta0 <- rep(0,ncol(x))
  if (missing(V1)) V1 <- diag(ncol(x))
  qtheta <- t(matrix(x,nrow=nrow(x))) - theta0
  qtheta <- qtheta %*% solve(n*g*V1) %*% t(qtheta)
  k <- -0.5*log(n*g*det(V1)) - lgamma(nu/2) + lgamma(ncol(x)/2) - .5*ncol(x)*log(pi)
  p1 <- ncol(x)
}
ans <- (k - .5*(nu+p1)*log(qtheta) -1/qtheta)
if (!logscale) { ans <- exp(ans); ans[is.na(ans)] <- 0 }
return(ans)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/dmom.r"
dmom <- function(x,V1=1,g=1,n=1,theta0,baseDensity='normal',nu=3) {

if (!(baseDensity %in% c('normal','t'))) stop("The only implemented baseDensity values are 'normal' and 't'")
if (missing(V1)) {
  if (is.vector(x)) V1 <- 1 else V1 <- diag(ncol(x))
}
if (missing(theta0)) {
  if (is.vector(V1)) theta0 <- 0 else theta0 <- rep(0,ncol(V1))
}
    
if (is.vector(V1)) {
  qtheta <- (x-theta0)^2/(n*g*V1)
  if (baseDensity=='normal') {
    ans <- qtheta*dnorm(x,theta0,sd=sqrt(n*g*V1))
  } else if (baseDensity=='t') {
    normct <- exp(lgamma(.5*(nu+1))-lgamma(.5*nu)-.5*log(nu*pi*n*g*V1))
    ans <- qtheta*normct*(1+qtheta/nu)^(-.5*(nu+1))*(nu-2)/nu
  }
} else {
  require(mvtnorm)
  qtheta <- mahalanobis(x,center=theta0,cov=n*g*V1)
  if (baseDensity=='normal') {
    ans <- qtheta*dmvnorm(x,mean=theta0,sigma=n*g*V1)/ncol(V1)
  } else if (baseDensity=='t') {
    ans <- qtheta*dmvt(x,delta=theta0,sigma=n*g*V1)*(nu-2)/(nu*length(theta0))
  }
}
return(ans)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/g2mode.r"
g2mode <- function(g,prior='iMom',nu=1,dim=1) {
  if (!(prior %in% c('normalMom','tMom','iMom'))) stop("Currently only prior=='normalMom', 'tMom' or 'iMom' are implemented")
  if (prior=='normalMom') {
    return(2*g)
  } else if (prior=='tMom') {
    return(g*2*nu/(nu-2+dim))
  } else if (prior=='iMom') {
    return(2*g/(nu+dim))
  }
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/imombf.lm.r"
imombf.lm <- function(lm1,coef,g,prior.mode,nu=1,theta0,method='adapt',nquant=100,B=10^5) {
if ((!missing(g)) & (!missing(prior.mode))) warning('Both g and prior.mode were specified. g will be ignored')
if ((missing(g)) & (missing(prior.mode))) stop('Either g or prior.mode must be specified')
if (missing(theta0)) theta0 <- rep(0,length(coef)) else if (length(theta0)!=length(coef)) stop('theta0 must have the same length as coef')
  
  thetahat <- coef(lm1)
  V <- summary(lm1)$cov.unscaled
  n <- length(lm1$residuals); p <- length(thetahat); p1 <- length(coef)
  if ((min(coef)<1) | (max(coef)>p)) stop('Non-valid value for coef. Use only values between 1 and the number of coefficients in lm1')
  ssr <- sum(residuals(lm1)^2); sr <- sqrt(ssr/(n-p))
  if (!missing(prior.mode)) g <- mode2g(prior.mode,prior='iMom')
  bf.imom <- imomunknown(thetahat[coef],V[coef,coef],n,nuisance.theta=p-p1,g=g,nu=nu,theta0=theta0,ssr=ssr,method=method,nquant=nquant,B=B)
  return(bf.imom)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/imombf.r"
imombf <- function(lm1,coef,g,prior.mode,nu=1,theta0,method='adapt',nquant=100,B=10^5) { UseMethod("imombf") }
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/imomknown.r"
imomknown <- function(theta1hat,V1,n,nuisance.theta,g=1,nu=1,theta0,sigma,method='adapt',B=10^5) {
if (missing(sigma)) stop('sigma must be specified')
if (missing(theta0)) theta0 <- rep(0,length(theta1hat))
f <- function(z) { ans <- (n*gi/z)^((nu+p1)/2) * exp(-n*gi/z); ans[z==0] <- 0; return(ans) }

p1 <- length(theta1hat)
l <- theta1hat-theta0; l <- matrix(l,nrow=1) %*% solve(V1) %*% matrix(l,ncol=1) / sigma^2 #noncentr param
m <- double(length(g))
if (method=='MC') {
  z <- rchisq(B,df=p1,ncp=l)
  for (i in 1:length(m)) { gi <- g[i]; m[i] <- mean(f(z)) }
} else if (method=='adapt') {
  f2 <- function(z) { return(f(z)*dchisq(z,df=p1,ncp=l)) }
  for (i in 1:length(m)) { gi <- g[i]; m[i] <- integrate(f2,0,Inf)$value }
} else {
  stop('method must be adapt or MC')
}
bf <- exp((p1/2)*log(2/(n*g)) + lgamma(p1/2)-lgamma(nu/2) + .5*l) * m
return(bf)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/imomunknown.r"
imomunknown <- function(theta1hat,V1,n,nuisance.theta,g=1,nu=1,theta0,ssr,method='adapt',nquant=100,B=10^5) {

fncp <- function(sigma2) {
  l <- theta1hat-theta0; l <- matrix(l,nrow=1) %*% solve(V1) %*% matrix(l,ncol=1) / sigma2
  l[l==Inf] <- exp(80)
  return(l)
}
f <- function(z) { ans <- (n*gi/z)^((nu+p1)/2) * exp(-n*gi/z); ans[z==0] <- 0; return(ans) }

if (missing(theta0)) theta0 <- rep(0,length(theta1hat))
p1 <- length(theta1hat)
m <- double(length(g))
if (method=='MC') {
  sigma2 <- 1/rgamma(B,(n-nuisance.theta)/2,ssr/2)
  l <- fncp(sigma2)
  z <- rchisq(B,df=p1,ncp=l)
  for (i in 1:length(m)) { gi <- g[i]; m[i] <- mean(f(z)) }
} else if (method=='adapt') {
  f2 <- function(z,z2) { return(f(z)*dchisq(z,df=p1,ncp=fncp(z2))) }
  qseq <- 1/qgamma((2*(1:nquant)-1)/(2*nquant),(n-nuisance.theta)/2,ssr/2)
  for (i in 1:length(m)) {
    m[i] <- 0; for (j in 1:nquant) { gi <- g[i]; m[i] <- m[i]+integrate(f2,0,Inf,z2=qseq[j])$value }
    m[i] <- m[i]/nquant
  }
} else {
  stop('method must be adapt or MC')
}
t1 <- theta1hat-theta0; t1 <- matrix(t1,nrow=1) %*% solve(V1) %*% matrix(t1,ncol=1) / ssr #noncentr param
bf <- exp((p1/2)*log(2/(n*g)) + lgamma(p1/2) - lgamma(nu/2) + ((n-nuisance.theta)/2)*log(1+t1)) * m
return(bf)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/mode2g.r"
mode2g <- function(prior.mode,prior='iMom',nu=1,dim=1) {
  if (!(prior %in% c('normalMom','tMom','iMom'))) stop("Currently only prior=='normalMom', 'tMom' or 'iMom' are implemented")
  if (prior=='normalMom') {
    return(prior.mode/2)
  } else if (prior=='tMom') {
    if (nu<3) stop('The tMom prior must have nu>2 degrees of freedom')
    return(prior.mode*(nu-2+dim)/(2*nu))    
  } else if (prior=='iMom') {
    return(prior.mode *(nu+dim)/2)
  }
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/mombf.lm.r"
mombf.lm <- function(lm1,coef,g,prior.mode,baseDensity='normal',nu=3,theta0,logbf=FALSE,B=10^5) {
 
if ((!missing(g)) & (!missing(prior.mode))) warning('Both g and prior.mode were specified. g will be ignored')
if ((missing(g)) & (missing(prior.mode))) stop('Either g or prior.mode must be specified')
if (missing(theta0)) theta0 <- rep(0,length(coef)) else if (length(theta0)!=length(coef)) stop('theta0 must have the same length as coef')
  
  thetahat <- coef(lm1)
  V <- summary(lm1)$cov.unscaled
  n <- length(lm1$residuals); p <- length(thetahat); p1 <- length(coef)
  if ((min(coef)<1) | (max(coef)>p)) stop('Non-valid value for coef. Use only values between 1 and the number of coefficients in lm1')
  ssr <- sum(residuals(lm1)^2); sr <- sqrt(ssr/(n-p))
  if (!missing(prior.mode)) g <- mode2g(prior.mode,prior=paste(baseDensity,'Mom',sep=''),nu=nu)
  if (baseDensity=='normal') {
    bf.mom <- momunknown(thetahat[coef],V[coef,coef],n=n,g=g,ssr=ssr,nuisance.theta=p-p1,theta0=theta0,logbf=logbf)
  } else if (baseDensity=='t') {
    randomg <- 1/rgamma(B,nu/2,nu/2)
    bf.mom <- momunknown(thetahat[coef],V[coef,coef],n=n,g=g*randomg,ssr=ssr,nuisance.theta=p-p1,theta0=theta0,logbf=logbf)
    bf.mom <- mean(bf.mom*g)*(nu-2)/nu
  }
  return(bf.mom)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/mombf.r"
mombf <- function(lm1,coef,g,prior.mode,baseDensity='normal',nu=3,theta0,logbf=FALSE,B=10^5) {
  if (!(baseDensity %in% c('normal','t'))) stop('baseDensity must either be normal or t')
  if (baseDensity=='t') {
    if (length(g)>1) stop("For baseDensity=='t' g is only allowed to have length 1")
    if (nu <= 2) stop('For t base density nu must be >=3, otherwise prior is improper')
    if (nu >= length(lm1$residuals)-length(coef)+2) warning(paste('Based on the amount of data available, nu should be <',length(lm1$residuals)-length(coef)+2,'in order to guarantee finite sample consistency'))
  }
  UseMethod("mombf")
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/momknown.r"
momknown <- function(theta1hat,V1,n,g=1,theta0,sigma,logbf=FALSE) {
if (missing(sigma)) stop('sigma must be specified')
if (missing(theta0)) theta0 <- rep(0,length(theta1hat))
p1 <- length(theta1hat)
l <- theta1hat-theta0; l <- matrix(l,nrow=1) %*% solve(V1) %*% matrix(l,ncol=1) * n*g/((1+n*g)*sigma^2) #noncentr param
muk <- p1+l
t1 <- matrix(theta1hat-theta0,nrow=1) %*% solve(V1) %*% matrix(theta1hat-theta0,ncol=1) * n*g/((1+n*g)*sigma^2)
bf <- .5*t1 + log(muk) - log(1+n*g) - (p1/2)*log(1+n*g) - log(p1)
if (!logbf) bf <- exp(bf)
return(bf)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/momunknown.r"
momunknown <- function(theta1hat,V1,n,nuisance.theta,g=1,theta0,ssr,logbf=FALSE) {
if (missing(theta0)) theta0 <- rep(0,length(theta1hat))
p1 <- length(theta1hat); p <- p1 + nuisance.theta
l <- theta1hat-theta0; l <- matrix(l,nrow=1) %*% solve(V1) %*% matrix(l,ncol=1)  
sigma2hat <- (ssr + l/(1+n*g))/(n-nuisance.theta)
muk <- p1+ l* n*g/((1+n*g)*sigma2hat)
bf <- (-(n-nuisance.theta)/2)*log(1+n*g*ssr/(ssr+l)) + log(muk) - log(1+n*g) + ((n-p)/2)*log(1+n*g) - log(p1)
if (!logbf) bf <- exp(bf)
return(bf)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/pimom.r"
pimom <- function(q,V1=1,g=1,n=1,nu=1) {

  z <- (.5*pgamma(n*g*V1/q^2,nu/2,1))
  return(z*(q<=0)+(1-z)*(q>0))
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/pmom.r"
pmom <- function(q,V1=1,g=1,n=1) {

  z <- .5-(pnorm(abs(q)/sqrt(n*V1*g)) - abs(q)/sqrt(2*pi*n*V1*g) * exp(-.5*q^2/(g*n*V1)) - .5)
  return(z*(q<=0)+(1-z)*(q>0))
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/priorp2g.r"
priorp2g <- function(priorp,q,nu=1,prior="iMom") {
  if (!(prior %in% c('normalMom','tMom','iMom'))) stop("Currently only prior=='normalMom', 'tMom' or 'iMom' are implemented")
  if (prior == "normalMom") {
    e <- function(logg) { return((1-2*pmom(-abs(q),g=exp(logg))-priorp[i])^2) }
    ans <- double(length(priorp))
    for (i in 1:length(priorp)) {
      ans[i] <- exp(nlminb(start=0,objective=e)$par)
    }
  } else if (prior == 'tMom') {
    stop("prior=='tMom' is not currently implemented in priorp2g")
  } else if (prior == "iMom") {
    p <- (1-priorp)/2
    ans <- qgamma(2*p,nu/2,1)*q^2
  }
  return(ans)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/qimom.r"
qimom <- function(p,V1=1,g=1,n=1,nu=1) {

  ans <- double(length(p))
  ans[p<=.5] <- -sqrt(n*g*V1/qgamma(2*p[p<=.5],nu/2,1))
  ans[p>.5] <- sqrt(n*g*V1/qgamma(2*(1-p[p>.5]),nu/2,1))
  return(ans)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/qmom.r"
qmom <- function(p,V1=1,g=1,n=1) {

  e <- function(q) { return((pmom(q,V1,g,n)-pneg)^2) }
  ans <- double(length(p))
  for (i in 1:length(p)) {
    pneg <- ifelse(p[i]<=.5,p[i],1-p[i])
    ans[i] <- nlminb(start=-1,objective=e)$par
    if (p[i]>.5) ans[i] <- -ans[i]
  }
  return(ans)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/zbfknown.r"
zbfknown <- function(theta1hat,V1,n,g=1,theta0,sigma,logbf=FALSE) {

if (missing(sigma)) stop('sigma must be specified')
if (missing(theta0)) theta0 <- rep(0,length(theta1hat))
p1 <- length(theta1hat)
l <- theta1hat-theta0; l <- matrix(l,nrow=1) %*% solve(V1) %*% matrix(l,ncol=1) * n*g/((1+n*g)*sigma^2) #noncentr param
muk <- p1+l
t1 <- matrix(theta1hat-theta0,nrow=1) %*% solve(V1) %*% matrix(theta1hat-theta0,ncol=1) * n*g/((1+n*g)*sigma^2)
bf <- .5*t1 - .5*p1*log(1+n*g)
if (!logbf) bf <- exp(bf)
return(bf)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/zbfunknown.r"
zbfunknown <- function(theta1hat,V1,n,nuisance.theta,g=1,theta0,ssr,logbf=FALSE) {

if (missing(theta0)) theta0 <- rep(0,length(theta1hat))
p1 <- length(theta1hat); p <- p1 + nuisance.theta
l <- theta1hat-theta0; l <- matrix(l,nrow=1) %*% solve(V1) %*% matrix(l,ncol=1)  
sigma2hat <- (ssr + l/(1+n*g))/(n-nuisance.theta)
muk <- p1+ l* n*g/((1+n*g)*sigma2hat)
bf <- (-(n-nuisance.theta)/2)*log(1+n*g*ssr/(ssr+l)) + ((n-p)/2)*log(1+n*g)
if (!logbf) bf <- exp(bf)
return(bf)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/zellnerbf.knownsig.r"
zellnerbf.knownsig <- function(theta1hat,V1,n,g=1,theta0,sigma,logbf=FALSE) {
#Bayes factor based on Zellner's g-prior for linear models (known sigma^2 case).
# - theta1hat: vector with estimated value of the coefficients that are to be tested
# - V1: submatrix of covariance corresponding to elements in theta1hat
# - n: sample size used to fit the model
# - g: prior parameter
# - theta0: hypothesized value for theta1hat (defaults to 0)
# - sigma: residual standard deviation, assumed to be known
# - logbf: if log==TRUE, the natural logarithm of the Bayes factor is returned
if (missing(sigma)) stop('sigma must be specified')
if (missing(theta0)) theta0 <- rep(0,length(theta1hat))
p1 <- length(theta1hat)
l <- theta1hat-theta0; l <- matrix(l,nrow=1) %*% solve(V1) %*% matrix(l,ncol=1) * n*g/((1+n*g)*sigma^2) #noncentr param
muk <- p1+l
t1 <- matrix(theta1hat-theta0,nrow=1) %*% solve(V1) %*% matrix(theta1hat-theta0,ncol=1) * n*g/((1+n*g)*sigma^2)
bf <- .5*t1
if (!logbf) bf <- exp(bf)
return(bf)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/zellnerbf.lm.r"
zellnerbf.lm <- function(lm1,coef,g,theta0,logbf=FALSE) {
  
if (missing(g)) stop('g must be specified')
if (missing(theta0)) theta0 <- rep(0,length(coef)) else if (length(theta0)!=length(coef)) stop('theta0 must have the same length as coef')
  
  thetahat <- coef(lm1)
  V <- summary(lm1)$cov.unscaled
  n <- length(lm1$residuals); p <- length(thetahat); p1 <- length(coef)
  if ((min(coef)<1) | (max(coef)>p)) stop('Non-valid value for coef. Use only values between 1 and the number of coefficients in lm1')
  ssr <- sum(residuals(lm1)^2); sr <- sqrt(ssr/(n-p))
  bf.zellner <- zbfunknown(thetahat[coef],V[coef,coef],n=n,nuisance.theta=p-p1,g=g,theta0=theta0,ssr=ssr,logbf=logbf)
  return(bf.zellner)
}
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/zellnerbf.r"
zellnerbf <- function(lm1,coef,g,theta0,logbf=FALSE) { UseMethod("zellnerbf") }
#line 1 "/private/tmp/RtmpQTS6su/R.INSTALL10d63af1/mombf/R/zellnerbf.unknownsig.r"
zellnerbf.unknownsig <- function(theta1hat,V1,n,nuisance.theta,g=1,theta0,ssr,logbf=FALSE) {
#Bayes factor based on Zellner's g-prior for linear models (unknown sigma^2 case). 
# - theta1hat: vector with estimated value of the coefficients that are to be tested
# - V1: submatrix of covariance corresponding to elements in theta1hat
# - n: sample size used to fit the model
# - g: prior parameter
# - theta0: hypothesized value for theta1hat (defaults to 0)
# - ssr: sum of squared residuals
# - logbf: if log==TRUE, the natural logarithm of the Bayes factor is returned
if (missing(theta0)) theta0 <- rep(0,length(theta1hat))
p1 <- length(theta1hat); p <- p1 + nuisance.theta
l <- theta1hat-theta0; l <- matrix(l,nrow=1) %*% solve(V1) %*% matrix(l,ncol=1)  
sigma2hat <- (ssr + l/(1+n*g))/(n-nuisance.theta)
muk <- p1+ l* n*g/((1+n*g)*sigma2hat)
bf <- (-(n-nuisance.theta)/2)*log(1+n*g*ssr/(ssr+l)) + ((n-p)/2)*log(1+n*g)
if (!logbf) bf <- exp(bf)
return(bf)
}
