Update logLik

haziqj · Jun 8, 2017 · a42007b · a42007b
1 parent 257136b
commit a42007b
Show file tree

Hide file tree

Showing 6 changed files with 53 additions and 53 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -2,7 +2,7 @@
 
 S3method(fitted,ipriorProbit)
 S3method(fitted,iprobitMult)
-S3method(logLik,ipriorProbit)
+S3method(logLik,iprobitMod)
 S3method(plot,ipriorProbit)
 S3method(plot,iprobitData)
 S3method(plot,iprobitMult)

diff --git a/R/EprodPhiZ.R b/R/EprodPhiZ.R
diff --git a/R/iprobit_helper.R b/R/iprobit_helper.R
@@ -36,3 +36,47 @@ iprobitSE <- function(y, eta, thing1 = NULL, thing0 = NULL) {
   var.ystar[y == 0] <- 1 - eta[y == 0] * thing0 + (thing0 ^ 2)
   sqrt(var.ystar)
 }
+
+# This is wrong
+# EprodPhiZ_1 <- function(mu, sigma) {
+#   Sigma <- outer(sigma, sigma, FUN = "*") + diag(1, length(mu))
+#   L.inv <- solve(t(chol(Sigma)))
+#   e <- L.inv %*% mu
+#   prod(pnorm(e))
+# }
+#
+# EprodPhiZ_2 <- function(mu, sigma) {
+#   Sigma <- outer(sigma, sigma, FUN = "*") + diag(1, length(mu))
+#   pmvnorm(upper = mu, sigma = Sigma)
+# }
+
+EprodPhiZ <- function(mu, sigma = rep(1, length(mu)), log = FALSE) {
+  res <- integrate(
+    function(z) {
+      tmp <- 0
+      for (i in seq_len(length(mu)))
+        tmp <- tmp + pnorm(z + mu[i], log.p = TRUE)
+      exp(tmp) * dnorm(z)
+    },
+    lower = -Inf, upper = Inf
+  )$value
+  ifelse(isTRUE(log), log(res), res)
+}
+
+# EprodPhiZ <- function(mu, sigma) {
+#   res1 <- EprodPhiZ_1(mu = mu, sigma = sigma)
+#   res2 <- EprodPhiZ_2(mu = mu, sigma = sigma)
+#   res3 <- EprodPhiZ_3(mu = mu, sigma = sigma)
+#
+#   list(res1, res2, res3)
+# }
+#
+# k <- 20
+# microbenchmark::microbenchmark(
+#   method1 = EprodPhiZ_1(rnorm(k), rep(1, k)),
+#   method2 = EprodPhiZ_2(rnorm(k), rep(1, k)),
+#   method3 = EprodPhiZ_3(rnorm(k), rep(1, k))
+# )
+#
+#
+# EprodPhiZ(rnorm(4), rep(1, 4))
diff --git a/R/logLik.R b/R/logLik.R
@@ -4,7 +4,7 @@
 #'
 #' @return The variational lower bound.
 #' @export
-logLik.ipriorProbit <- function(x) {
+logLik.iprobitMod <- function(x) {
   lb <- x$lower.bound[!is.na(x$lower.bound)]
   lb <- lb[length(lb)]
   class(lb) <- "iprobitLowerBound"
@@ -14,4 +14,4 @@ logLik.ipriorProbit <- function(x) {
 #' @export
 print.iprobitLowerBound <- function(x) {
   cat("Lower bound =", x)
-}
+}
diff --git a/R/iprobit_print_and_summary.R → R/print_and_summary.R b/R/iprobit_print_and_summary.R → R/print_and_summary.R
@@ -1,12 +1,12 @@
 interceptAndLambdaNames <- function(x) {
   Intercept <- x$alpha
   if (length(Intercept) > 1)
-    names(Intercept) <- paste0("Intercept[", seq_along(Intercept), "] ")
+    names(Intercept) <- paste0("alpha[", seq_along(Intercept), "]")
   else
-    names(Intercept) <- "Intercept"
+    names(Intercept) <- "alpha"
   lambda <- x$lambda
   if (length(lambda) > 1)
-    names(lambda) <- paste0("lambda[", seq_along(lambda), "] ")
+    names(lambda) <- paste0("lambda[", seq_along(lambda), "]")
   else
     names(lambda) <- "lambda"
   c(Intercept, lambda)

diff --git a/man/logLik.ipriorProbit.Rd → man/logLik.iprobitMod.Rd b/man/logLik.ipriorProbit.Rd → man/logLik.iprobitMod.Rd