[r-cran-mcmcpack] 22/90: Imported Upstream version 0.8-2
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:35 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-mcmcpack.
commit 19673ae3f8c0c9f95f6d81f9466f25476885d425
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:09 2016 +0100
Imported Upstream version 0.8-2
---
DESCRIPTION | 6 +-
HISTORY | 8 +
R/MCMClogit.R | 28 +--
R/MCMCmetrop1R.R | 2 +-
R/MCMCpoisson.R | 25 ++-
R/MCMCpoissonChangepoint.R | 420 ++++++++++++++++-----------------------------
R/MCMCprobit.R | 24 ++-
R/MCMCregress.R | 35 ++--
R/hidden.R | 19 ++
9 files changed, 248 insertions(+), 319 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 5a1e4b2..de09038 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: MCMCpack
-Version: 0.8-1
-Date: 2007-1-11
+Version: 0.8-2
+Date: 2007-4-23
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin <admartin at wustl.edu>, and
Kevin M. Quinn <kevin_quinn at harvard.edu>
@@ -17,4 +17,4 @@ Description: This package contains functions to perform Bayesian
sampling algorithm, and tools for visualization.
License: GPL version 2
URL: http://mcmcpack.wustl.edu
-Packaged: Thu Jan 11 17:48:15 2007; adm
+Packaged: Mon Apr 23 15:23:25 2007; adm
diff --git a/HISTORY b/HISTORY
index 53637bd..26e2d0b 100644
--- a/HISTORY
+++ b/HISTORY
@@ -2,6 +2,14 @@
// Changes and Bug Fixes
//
+0.8-1 to 0.8-2
+ * models with multivariate normal priors now check for symmetry and positive
+ semi-definiteness of precision matrix
+ * fixed bug in how models with marginal likelihood calculations check for
+ prior propriety [thanks to Gary Rosner for spotting this]
+ * fixed bug in how MCMCmetrop1R handled a singular Hessian
+ [thanks to Piers Dunstan for spotting this]
+
0.7-4 to 0.8-1
* added MCMCPoissonChangepoint() model (authored by Jong Hee Park)
* added two plot methods for changepoint models:
diff --git a/R/MCMClogit.R b/R/MCMClogit.R
index a30cb4c..1d09c8d 100644
--- a/R/MCMClogit.R
+++ b/R/MCMClogit.R
@@ -22,17 +22,7 @@
cl <- match.call()
- ## get marginal likelihood argument
- marginal.likelihood <- match.arg(marginal.likelihood)
- if (max(B0==0) == 1 & is.null(user.prior.density)){
- if (marginal.likelihood != "none"){
- warning("Cannot calculate marginal likelihood with improper prior\n")
- marginal.likelihood <- "none"
- }
- }
- logmarglike <- NULL
-
-
+
## seeds
seeds <- form.seeds(seed)
lecuyer <- seeds[[1]]
@@ -52,6 +42,22 @@
b0 <- mvn.prior[[1]]
B0 <- mvn.prior[[2]]
+
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+ B0.eigenvalues <- eigen(B0)$values
+ if (min(B0.eigenvalues) < 0){
+ stop("B0 is not positive semi-definite.\nPlease respecify and call again.\n")
+ }
+ if (isTRUE(all.equal(min(B0.eigenvalues), 0))){
+ if (marginal.likelihood != "none"){
+ warning("Cannot calculate marginal likelihood with improper prior\n")
+ marginal.likelihood <- "none"
+ }
+ }
+ logmarglike <- NULL
+
+
## setup the environment so that fun can see the things passed as ...
userfun <- function(ttt) user.prior.density(ttt, ...)
my.env <- environment(fun = userfun)
diff --git a/R/MCMCmetrop1R.R b/R/MCMCmetrop1R.R
index 267f10b..41ba0f3 100644
--- a/R/MCMCmetrop1R.R
+++ b/R/MCMCmetrop1R.R
@@ -70,7 +70,7 @@
hess.new <- opt.out$hessian
hess.flag <- 0
if (force.samp==TRUE){
- if (max(diag(optim.out$hessian)==0)){
+ if (max(diag(opt.out$hessian)==0)){
for (i in 1:nrow(hess.new)){
if (hess.new[i,i] == 0){
hess.new[i,i] <- -1e-6
diff --git a/R/MCMCpoisson.R b/R/MCMCpoisson.R
index 16253e7..e86b9fa 100644
--- a/R/MCMCpoisson.R
+++ b/R/MCMCpoisson.R
@@ -19,15 +19,6 @@
cl <- match.call()
- ## get marginal likelihood argument
- marginal.likelihood <- match.arg(marginal.likelihood)
- if (max(B0==0) == 1){
- if (marginal.likelihood != "none"){
- warning("Cannot calculate marginal likelihood with improper prior\n")
- marginal.likelihood <- "none"
- }
- }
- logmarglike <- NULL
## seeds
seeds <- form.seeds(seed)
@@ -47,6 +38,22 @@
mvn.prior <- form.mvn.prior(b0, B0, K)
b0 <- mvn.prior[[1]]
B0 <- mvn.prior[[2]]
+
+
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+ B0.eigenvalues <- eigen(B0)$values
+ if (min(B0.eigenvalues) < 0){
+ stop("B0 is not positive semi-definite.\nPlease respecify and call again.\n")
+ }
+ if (isTRUE(all.equal(min(B0.eigenvalues), 0))){
+ if (marginal.likelihood != "none"){
+ warning("Cannot calculate marginal likelihood with improper prior\n")
+ marginal.likelihood <- "none"
+ }
+ }
+ logmarglike <- NULL
+
## form the tuning parameter
tune <- vector.tune(tune, K)
diff --git a/R/MCMCpoissonChangepoint.R b/R/MCMCpoissonChangepoint.R
index d4e9579..4669cb7 100644
--- a/R/MCMCpoissonChangepoint.R
+++ b/R/MCMCpoissonChangepoint.R
@@ -1,272 +1,148 @@
-################################
-## Poisson Changepoint Model ##
-################################
-
-## This version does not allow covariates
-
-## Noninformative prior is not acceptable
-## in case of marginal likelihood computation
-
-## Originally written 12/23/2005 Jong Hee Park
-## Modified 02/13/2006 Jong Hee Park
-## Modified 10/21/2006 Jong Hee Park
-## Modified 12/13/2006 Jong Hee Park for JStatSoft article
-
-## some utility functions added
-## Jong Hee Park 07/14/2006
-
-## two plot functions added
-## Jong Hee Park 12/13/2006
-
-##############################################################
-## Andrew and Kevin, followings are a list of things to do
-##############################################################
-
-## 1. When you write a C code for MCMCpoissonChangepoint,
-## you need to export "prob.state" in addition to mcmc output.
-## This "prob.state" contains posterior probabilities of each state
-## and is essential to draw plots and posterior inference.
-
-## 2. I provide two special plot functions for changepoint models:
-## "plot.post.state" and "plot.post.changepoint."
-## These plot functions need to be changed based on the changes in
-## MCMCpoissonChangepoint() output.
-
-## 3. We need to change outputs of MCMCpoissonChangepoint make it accessible
-## by BayesFactor().
-
-## 4. All helper functions should be straightforward. "rdirichlet.cp" and
-## "ddirichlet.cp" might be redundant because they're simple a beta
-## distribution. In the future, I will replace these into beta functions
-## but for the time being, use these functions. Assigning zeros in the right
-## place is trickier than I thought. (I should have written these using
-## beta density in the first place...)
-
-
-
-#########################################################
-"MCMCpoissonChangepoint"<-
-#########################################################
- function(data, m = 1, burnin = 1000, mcmc = 1000, thin = 1,
- verbose = 0, seed = NA, c0, d0, a = NULL, b = NULL,
- marginal.likelihood = c("none", "Chib95"), ...) {
-
- ####################################################
- ## Arguments
- ## m : the number of changepoint
- ## c0 and d0: gamma prior for lambda. NO DEFAULT values
- ## a and b: beta prior for transition probabilities
- ## By default, the expected duration is computed and
- ## corresponding a and b values are assigned. The expected
- ## duration is the sample period divided by the number of states.
- ## ex) 300 years and 2 changepoints (3 states)
- ## expected.duration <- 300/(2+1)
- ## b <- 0.1
- ## a <- b*expected.duration
- ####################################################
-
- # check iteration parameters
- check.mcmc.parameters(burnin, mcmc, thin)
- totiter <- mcmc + burnin
- cl <- match.call()
-
- ## ## seeds
- ## seeds <- form.seeds(seed)
- ## lecuyer <- seeds[[1]]
- ## seed.array <- seeds[[2]]
- ## lecuyer.stream <- seeds[[3]]
- if(!is.na(seed)) set.seed(seed)
-
- ## sample size
- y <- data
- n <- length(y)
-
- ## prior
- A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
-
- ## get marginal likelihood argument
- marginal.likelihood <- match.arg(marginal.likelihood)
-
- ## storage matrices
- lambda.store <- matrix(NA, mcmc/thin, m+1)
- P.store <- matrix(NA, mcmc/thin, (m+1)^2)
- ps.store <- matrix(0, n, m+1)
- s1.store <- matrix(NA, mcmc/thin, n)
- py <- rep(0, m+1)
- pdf.P.store <- matrix(NA, mcmc, m+1)
- lambda1 <- rep(NA, m+1)
- P1 <- matrix(NA, m+1, m+1)
-
- ## starting values
- lambda0 <- runif(m+1)
- P0 <- trans.mat.prior(m=m, n=n, a=0.9, b=0.1)
-
- ####################################################################
- ## MCMC iteration starts!
- ####################################################################
-
- for (iter in 1:totiter){
-
- ####################################################################
- ## Step 1: Sampling S
- ####################################################################
- state.out <- Poisson.state.sampler(m=m, y=y, lambda=lambda0, P=P0)
- s1 <- state.out$s1
- ps1<- state.out$ps1
-
- ####################################################################
- ## Step 2: Sampling lambda
- ####################################################################
- for (j in 1:(m+1)){
- ej <- as.numeric(s1==j)
- yj <- y[ej==1]
- nj <- length(yj)
- c1 <- sum(yj) + c0
- d1 <- nj + d0
- lambda1[j] <- rgamma(1, c1, d1)
- }
-
- ####################################################################
- ## Step 3: Sampling P
- ####################################################################
- switch <- switchg(s1)
- for (j in 1:(m+1)){
- switch1 <- A0[j,] + switch[j,]
- pj <- rdirichlet.cp(1, switch1)
- P1[j,] <- pj
- }
-
- ## update
- lambda0 <- lambda1
- P0 <- P1
-
- ####################################################################
- ## end of one iteration
- ####################################################################
-
- ## store
- if (iter > burnin && (iter %% thin == 0)) {
- lambda.store[iter-burnin,] <- lambda1
- P.store[iter-burnin,] <- as.vector(t(P1))
- s1.store[iter-burnin,] <- s1
- ps.store <- ps.store + ps1
- }
-
- ## report
- if(verbose > 0 && iter%%verbose == 0){
- cat("----------------------------------------------",'\n')
- cat("iteration = ", iter, '\n')
- cat("lambda = ", lambda1, '\n')
- cat("Transition Matrix", '\n')
- for(i in 1:(m+1))
- cat(paste("", P1[i,]), fill=TRUE, labels=paste("{",i,"}:", sep=""), sep=",")
- }
-
- } ## end of MCMC sampling
-
- ## marginal likelihood calculation if Chib == TRUE
- if (marginal.likelihood == "Chib95"){
-
- ############################################
- ## Bayes Factor Calculation Starts ##
- ############################################
- lambda.st <- apply(lambda.store, 2, mean)
- P.vec.st <- apply(P.store, 2, mean)
- P.st <- t(matrix(P.vec.st, m+1, m+1))
-
- #################################################################
- ## 1. pdf.lambda
- #################################################################
- density.lambda <- matrix(NA, mcmc, m+1)
- for (i in 1:mcmc){
- for (j in 1:(m+1)){
- # compute proposal density
- ej <- as.numeric(s1.store[i,]==j)
- yj <- y[ej==1]
- nj <- length(yj)
- c1 <- sum(yj) + c0
- d1 <- nj + d0
- density.lambda[i, j] <- dgamma(lambda.st[j], c1, d1)
- }
- }
- pdf.lambda <- log(prod(apply(density.lambda, 2, mean)))
-
- ######################################################################
- ## 2. pdf.P
- ######################################################################
- for (g in 1:mcmc){
- state.out <- Poisson.state.sampler(m=m, y=y, lambda=lambda.st, P=P0)
- s1 <- state.out$s1
- ps1 <- state.out$ps1
- switch <- switchg(s1)
-
- for (j in 1:(m+1)){
- switch1 <- A0[j,] + switch[j,]
- pj <- rdirichlet.cp(1, switch1)
- P1[j,] <- pj
- pdf.P.store[g,j]<- ddirichlet.cp(P.st[j,], switch1)
- }
- P0 <- P1
- }
- pdf.P <- log(prod(apply(pdf.P.store, 2, mean)))
-
- ####################################################################
- ## likelihood
- ####################################################################
- F <- matrix(NA, n, m+1)
- like<- rep(NA, n)
- pr1 <- c(1,rep(0,m))
- for (t in 1:n){
- py <- sapply(c(1:(m+1)), function(i){poisson.pdf(y[t], lambda.st[i])})
- if(t==1) {pstyt1 = pr1}
- else {pstyt1 <- F[t-1,]%*%P.st}
- unnorm.pstyt <- pstyt1*py
- pstyt <- unnorm.pstyt/sum(unnorm.pstyt)
- F[t,] <- pstyt
- like[t] <- sum(unnorm.pstyt)
- }
- loglik <- sum(log(like))
-
- ####################################################################
- ## prior ordinates
- ####################################################################
- nprior.lambda<- nprior.P <- rep(NA, m+1)
- nprior.lambda<- sapply(c(1:(m+1)), function(i){dgamma(lambda.st[i], c0, d0, log=TRUE)})
- nprior.P <- sapply(c(1:(m+1)), function(i){log(ddirichlet.cp(P.st[i,], A0[i,]))})
-
- prior.lambda <- sum(nprior.lambda)
- prior.P <- sum(nprior.P)
-
- ############################################
- ## Marginal Likelihood ##
- ############################################
- numerator <- loglik + prior.lambda + prior.P
- denominator <- pdf.lambda + pdf.P
- logmarglike <- numerator - denominator
-
- ## print
- if(verbose > 0){
- cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
- cat("Log Marginal Likelihood\n")
- cat("-------------------------------------------------",'\n')
- cat("log(marglike)= ", logmarglike, '\n')
- cat("log(likelihood)= ", loglik, '\n')
- cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
- }
- } ## end of marginal likelihood calculation
-
- else {marginal <- loglik <- NULL}
-
- ## return output
- output <- as.mcmc(lambda.store)
- varnames(output) <- paste("lambda.", 1:(m+1), sep = "")
- attr(output, "title") <- "MCMCpoissonChangepoint Posterior Sample"
- attr(output, "y") <- data
- attr(output, "call") <- cl
- attr(output, "logmarglike") <- logmarglike
- attr(output, "loglik") <- loglik
- attr(output, "prob.state") <- ps.store/mcmc
- return(output)
-
- }## end of MCMC function
-
+"MCMCpoissonChangepoint" <-
+ function (data, m = 1, burnin = 1000, mcmc = 1000, thin = 1, verbose = 0,
+ seed = NA, c0, d0, a = NULL, b = NULL, marginal.likelihood = c("none", "Chib95"), ...)
+{
+ check.mcmc.parameters(burnin, mcmc, thin)
+ totiter <- mcmc + burnin
+ cl <- match.call()
+ if (!is.na(seed))
+ set.seed(seed)
+ y <- data
+ n <- length(y)
+ A0 <- trans.mat.prior(m = m, n = n, a = a, b = b)
+ marginal.likelihood <- match.arg(marginal.likelihood)
+ lambda.store <- matrix(NA, mcmc/thin, m + 1)
+ P.store <- matrix(NA, mcmc/thin, (m + 1)^2)
+ ps.store <- matrix(0, n, m + 1)
+ s1.store <- matrix(NA, mcmc/thin, n)
+ py <- rep(0, m + 1)
+ pdf.P.store <- matrix(NA, mcmc/thin, m + 1)
+ lambda1 <- rep(NA, m + 1)
+ P1 <- matrix(NA, m + 1, m + 1)
+ lambda0 <- runif(m + 1)
+ P0 <- trans.mat.prior(m = m, n = n, a = 0.9, b = 0.1)
+ for (iter in 1:totiter) {
+ state.out <- Poisson.state.sampler(m = m, y = y, lambda = lambda0,
+ P = P0)
+ s1 <- state.out$s1
+ ps1 <- state.out$ps1
+ for (j in 1:(m + 1)) {
+ ej <- as.numeric(s1 == j)
+ yj <- y[ej == 1]
+ nj <- length(yj)
+ c1 <- sum(yj) + c0
+ d1 <- nj + d0
+ lambda1[j] <- rgamma(1, c1, d1)
+ }
+ switch <- switchg(s1)
+ for (j in 1:(m + 1)) {
+ switch1 <- A0[j, ] + switch[j, ]
+ pj <- rdirichlet.cp(1, switch1)
+ P1[j, ] <- pj
+ }
+ lambda0 <- lambda1
+ P0 <- P1
+ if (iter > burnin && (iter%%thin == 0)) {
+ lambda.store[(iter - burnin)/thin, ] <- lambda1
+ P.store[(iter - burnin)/thin, ] <- as.vector(t(P1))
+ s1.store[(iter - burnin)/thin, ] <- s1
+ ps.store <- ps.store + ps1
+ }
+ if (verbose > 0 && iter%%verbose == 0) {
+ cat("----------------------------------------------",
+ "\n")
+ cat("iteration = ", iter, "\n")
+ cat("lambda = ", lambda1, "\n")
+ cat("Transition Matrix", "\n")
+ for (i in 1:(m + 1)) cat(paste("", P1[i, ]), fill = TRUE,
+ labels = paste("{", i, "}:", sep = ""), sep = ",")
+ }
+ }
+ if (marginal.likelihood == "Chib95") {
+ lambda.st <- apply(lambda.store, 2, mean)
+ P.vec.st <- apply(P.store, 2, mean)
+ P.st <- t(matrix(P.vec.st, m + 1, m + 1))
+ density.lambda <- matrix(NA, (mcmc/thin), m + 1)
+ for (i in 1:(mcmc/thin)) {
+ for (j in 1:(m + 1)) {
+ ej <- as.numeric(s1.store[i, ] == j)
+ yj <- y[ej == 1]
+ nj <- length(yj)
+ c1 <- sum(yj) + c0
+ d1 <- nj + d0
+ density.lambda[i, j] <- dgamma(lambda.st[j],
+ c1, d1)
+ }
+ }
+ pdf.lambda <- log(prod(apply(density.lambda, 2, mean)))
+ for (g in 1:(mcmc/thin)) {
+ state.out <- Poisson.state.sampler(m = m, y = y,
+ lambda = lambda.st, P = P0)
+ s1 <- state.out$s1
+ ps1 <- state.out$ps1
+ switch <- switchg(s1)
+ for (j in 1:(m + 1)) {
+ switch1 <- A0[j, ] + switch[j, ]
+ pj <- rdirichlet.cp(1, switch1)
+ P1[j, ] <- pj
+ pdf.P.store[g, j] <- ddirichlet.cp(P.st[j, ],
+ switch1)
+ }
+ P0 <- P1
+ }
+ pdf.P <- log(prod(apply(pdf.P.store, 2, mean)))
+ F <- matrix(NA, n, m + 1)
+ like <- rep(NA, n)
+ pr1 <- c(1, rep(0, m))
+ for (t in 1:n) {
+ py <- sapply(c(1:(m + 1)), function(i) {
+ poisson.pdf(y[t], lambda.st[i])
+ })
+ if (t == 1) {
+ pstyt1 = pr1
+ }
+ else {
+ pstyt1 <- F[t - 1, ] %*% P.st
+ }
+ unnorm.pstyt <- pstyt1 * py
+ pstyt <- unnorm.pstyt/sum(unnorm.pstyt)
+ F[t, ] <- pstyt
+ like[t] <- sum(unnorm.pstyt)
+ }
+ loglik <- sum(log(like))
+ nprior.lambda <- nprior.P <- rep(NA, m + 1)
+ nprior.lambda <- sapply(c(1:(m + 1)), function(i) {
+ dgamma(lambda.st[i], c0, d0, log = TRUE)
+ })
+ nprior.P <- sapply(c(1:(m + 1)), function(i) {
+ log(ddirichlet.cp(P.st[i, ], A0[i, ]))
+ })
+ prior.lambda <- sum(nprior.lambda)
+ prior.P <- sum(nprior.P)
+ numerator <- loglik + prior.lambda + prior.P
+ denominator <- pdf.lambda + pdf.P
+ logmarglike <- numerator - denominator
+ if (verbose > 0) {
+ cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
+ cat("Log Marginal Likelihood\n")
+ cat("-------------------------------------------------",
+ "\n")
+ cat("log(marglike)= ", logmarglike, "\n")
+ cat("log(likelihood)= ", loglik, "\n")
+ cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
+ }
+ }
+ else {
+ logmarglike <- marginal <- loglik <- NULL
+ }
+ output <- as.mcmc(lambda.store)
+ varnames(output) <- paste("lambda.", 1:(m + 1), sep = "")
+ attr(output, "title") <- "MCMCpoissonChangepoint Posterior Sample"
+ attr(output, "y") <- data
+ attr(output, "call") <- cl
+ attr(output, "logmarglike") <- logmarglike
+ attr(output, "loglik") <- loglik
+ attr(output, "prob.state") <- ps.store/(mcmc/thin)
+ return(output)
+}
+
diff --git a/R/MCMCprobit.R b/R/MCMCprobit.R
index dde2e9a..395b47c 100644
--- a/R/MCMCprobit.R
+++ b/R/MCMCprobit.R
@@ -19,15 +19,6 @@
cl <- match.call()
- ## get marginal likelihood argument
- marginal.likelihood <- match.arg(marginal.likelihood)
- if (max(B0==0) == 1){
- if (marginal.likelihood != "none"){
- warning("Cannot calculate marginal likelihood with improper prior\n")
- marginal.likelihood <- "none"
- }
- }
- logmarglike <- NULL
## seeds
seeds <- form.seeds(seed)
@@ -49,6 +40,21 @@
b0 <- mvn.prior[[1]]
B0 <- mvn.prior[[2]]
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+ B0.eigenvalues <- eigen(B0)$values
+ if (min(B0.eigenvalues) < 0){
+ stop("B0 is not positive semi-definite.\nPlease respecify and call again.\n")
+ }
+ if (isTRUE(all.equal(min(B0.eigenvalues), 0))){
+ if (marginal.likelihood != "none"){
+ warning("Cannot calculate marginal likelihood with improper prior\n")
+ marginal.likelihood <- "none"
+ }
+ }
+ logmarglike <- NULL
+
+
## residuals setup
resvec <- NULL
if (is.logical(bayes.resid) && bayes.resid==TRUE){
diff --git a/R/MCMCregress.R b/R/MCMCregress.R
index 511d435..fd92ef5 100644
--- a/R/MCMCregress.R
+++ b/R/MCMCregress.R
@@ -21,19 +21,6 @@
cl <- match.call()
- ## get marginal likelihood argument
- marginal.likelihood <- match.arg(marginal.likelihood)
- if (max(B0==0) == 1){
- if (marginal.likelihood != "none"){
- warning("Cannot calculate marginal likelihood with improper prior\n")
- marginal.likelihood <- "none"
- }
- }
- logmarglike <- NULL
- chib <- 0
- if (marginal.likelihood == "Chib95"){
- chib <- 1
- }
## seeds
seeds <- form.seeds(seed)
@@ -54,7 +41,27 @@
b0 <- mvn.prior[[1]]
B0 <- mvn.prior[[2]]
check.ig.prior(c0, d0)
-
+
+
+ ## get marginal likelihood argument
+ marginal.likelihood <- match.arg(marginal.likelihood)
+ B0.eigenvalues <- eigen(B0)$values
+ if (min(B0.eigenvalues) < 0){
+ stop("B0 is not positive semi-definite.\nPlease respecify and call again.\n")
+ }
+ if (isTRUE(all.equal(min(B0.eigenvalues), 0))){
+ if (marginal.likelihood != "none"){
+ warning("Cannot calculate marginal likelihood with improper prior\n")
+ marginal.likelihood <- "none"
+ }
+ }
+ logmarglike <- NULL
+ chib <- 0
+ if (marginal.likelihood == "Chib95"){
+ chib <- 1
+ }
+
+
## define holder for posterior sample
sample <- matrix(data=0, mcmc/thin, K+1)
diff --git a/R/hidden.R b/R/hidden.R
index 7543a70..c481750 100644
--- a/R/hidden.R
+++ b/R/hidden.R
@@ -593,6 +593,9 @@
# prior precision
if(is.null(dim(B0))) {
+ if (length(B0) > K){
+ stop("B0 was passed as a vector longer than K.\nB0 must be either a scalar or a matrix.\nPlease respecify and call ", calling.function(), " again.\n", call.=FALSE)
+ }
B0 <- B0 * diag(K)
}
if((dim(B0)[1] != K) || (dim(B0)[2] != K)) {
@@ -600,6 +603,22 @@
stop("Please respecify and call ", calling.function(), " again.\n",
call.=FALSE)
}
+
+ ## check B0 for symmetry
+ symproblem <- FALSE
+ for (i in 1:K){
+ for (j in i:K){
+ if (B0[i,j] != B0[j,i]){
+ symproblem <- TRUE
+ }
+ }
+ }
+ if (symproblem){
+ cat("B0 is not symmetric.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+
return(list(b0,B0))
}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-mcmcpack.git
More information about the debian-science-commits
mailing list