diff options
author | Ken Kellner <ken@kenkellner.com> | 2023-12-06 21:01:42 -0500 |
---|---|---|
committer | Ken Kellner <ken@kenkellner.com> | 2023-12-06 21:04:06 -0500 |
commit | 3628ba7cadad4c09a8b306e37267bc6671f7636d (patch) | |
tree | f66cde06d18bd38d3602d83f626f421b5d8e7239 | |
parent | 030c81cf866c40b7bfa9bafbfbb1d9f2c44d73c5 (diff) |
Fix #30 and add plot and jags.View tests
-rw-r--r-- | R/densityplot.R | 55 | ||||
-rw-r--r-- | R/jags_View.R | 2 | ||||
-rw-r--r-- | R/traceplot.R | 12 | ||||
-rw-r--r-- | R/whiskerplot.R | 10 | ||||
-rw-r--r-- | inst/tinytest/test_jags.R | 14 | ||||
-rw-r--r-- | inst/tinytest/test_plots.R | 49 |
6 files changed, 107 insertions, 35 deletions
diff --git a/R/densityplot.R b/R/densityplot.R index b350f0a..2f7209b 100644 --- a/R/densityplot.R +++ b/R/densityplot.R @@ -24,36 +24,41 @@ param_density <- function(x, parameter, m_labels=FALSE){ #Get samples vals <- mcmc_to_mat(x$samples[,parameter]) - - # Get bandwidth, one value for all chains - bw <- mean(apply(vals, 2, stats::bw.nrd0)) + + if(any(is.na(vals))){ + plot(1:nrow(vals), rep(0, nrow(vals)), type='n', + xlab='Value', ylab='Density', main=paste('Density of',parameter)) + } else { + # Get bandwidth, one value for all chains + bw <- mean(apply(vals, 2, stats::bw.nrd0)) - from <- min(vals) - 3*bw # these are 'density's defaults... - to <- max(vals) + 3*bw # ... use these if no constraints - xx <- vals - mult <- 1 + from <- min(vals) - 3*bw # these are 'density's defaults... + to <- max(vals) + 3*bw # ... use these if no constraints + xx <- vals + mult <- 1 - # Check for non-negative constraint or probability - if (min(vals) >= 0 && min(vals) < 2 * bw) { # it's non-negative - from <- 0 - xx <- rbind(vals, -vals) - mult <- 2 - } - if (min(vals) >= 0 && max(vals) <= 1 && + # Check for non-negative constraint or probability + if (min(vals) >= 0 && min(vals) < 2 * bw) { # it's non-negative + from <- 0 + xx <- rbind(vals, -vals) + mult <- 2 + } + if (min(vals) >= 0 && max(vals) <= 1 && (min(vals) < 2 * bw || 1 - max(vals) < 2 * bw)) { # it's a probability - xx <- rbind(vals, -vals, 2-vals) - mult <- 3 - to <- 1 - } + xx <- rbind(vals, -vals, 2-vals) + mult <- 3 + to <- 1 + } - # Get densities - dens <- apply(xx, 2, function(x) stats::density(x, bw=bw, from=from, to=to)$y) * mult + # Get densities + dens <- apply(xx, 2, function(x) stats::density(x, bw=bw, from=from, to=to)$y) * mult - # Draw plot - x <- seq(from, to, length=nrow(dens)) - cols <- grDevices::rainbow(ncol(vals)) - graphics::matplot(x, dens, type='l', lty=1, col=cols, - xlab='Value', ylab='Density', main=paste('Density of',parameter)) + # Draw plot + x <- seq(from, to, length=nrow(dens)) + cols <- grDevices::rainbow(ncol(vals)) + graphics::matplot(x, dens, type='l', lty=1, col=cols, + xlab='Value', ylab='Density', main=paste('Density of',parameter)) + } #Add margin labels if necessary if(m_labels){ diff --git a/R/jags_View.R b/R/jags_View.R index 0bae3e0..4ab0460 100644 --- a/R/jags_View.R +++ b/R/jags_View.R @@ -8,7 +8,7 @@ jags.View <- function(x, title, digits=3){ } # Organize columns if(x$mcmc.info$n.chains!=1){y = x$summary[,c(1,2,3,5,7,10,11,8,9)] - } else {y = x$summary[,c(1,2,3,5,7,8,9)]} + } else {y = x$summary[,c(1,2,3,5,7,10,11)]} z <- as.data.frame(round(as.matrix(y),digits)) if(is.vector(y)){ z <- as.data.frame(t(z)) diff --git a/R/traceplot.R b/R/traceplot.R index 755d3aa..66069fb 100644 --- a/R/traceplot.R +++ b/R/traceplot.R @@ -30,9 +30,16 @@ param_trace <- function(x, parameter, m_labels=FALSE){ #Draw plot cols <- grDevices::rainbow(ncol(vals)) - graphics::matplot(1:nrow(vals), vals, type='l', lty=1, col=cols, + + if(all(is.na(vals))){ + plot(1:nrow(vals), rep(0, nrow(vals)), type='n', + xlab="Iterations", ylab="Value", + main=bquote(.(parameter)*","~hat(R) == .(Rhat))) + } else { + graphics::matplot(1:nrow(vals), vals, type='l', lty=1, col=cols, xlab='Iterations', ylab='Value', main=bquote(.(parameter)*","~hat(R) == .(Rhat))) + } #Add margin labels if necessary if(m_labels){ @@ -40,6 +47,3 @@ param_trace <- function(x, parameter, m_labels=FALSE){ graphics::mtext("Value", side=2, line=1.5, outer=TRUE) } } - -#General function for setting up plots -# get_plot_info now has its own file diff --git a/R/whiskerplot.R b/R/whiskerplot.R index 8b09bf1..0e85236 100644 --- a/R/whiskerplot.R +++ b/R/whiskerplot.R @@ -27,11 +27,15 @@ whiskerplot <- function(x,parameters,quantiles=c(0.025,0.975), #Plot parameter means n <- length(parameters) - graphics::plot(1:n, post_stats[1,], xaxt="n", - ylim=range(post_stats), xlim=c(0,n+1), + graphics::plot(1:n, post_stats[1,], xaxt="n", type='n', + ylim=range(post_stats, na.rm=TRUE), xlim=c(0,n+1), xlab="Parameters", ylab=paste0('Parameter mean and quantiles (',quantiles[1], - ' - ',quantiles[2],')'), pch=19, cex=1.5, ...) + ' - ',quantiles[2],')'), ...) + for (i in 1:n){ + if(any(is.na(post_stats[,i]))) next + graphics::points(i, post_stats[1,i], cex=1.5, pch=19) + } graphics::axis(side=1, at=1:n, labels=parameters) graphics::box() diff --git a/inst/tinytest/test_jags.R b/inst/tinytest/test_jags.R index 2218dfc..be64ffb 100644 --- a/inst/tinytest/test_jags.R +++ b/inst/tinytest/test_jags.R @@ -193,6 +193,18 @@ expect_true(all(is.na(out$summary[,"Rhat"]))) expect_true(all(is.na(out$summary[,"n.eff"]))) expect_true(all(out$summary["alpha",3:7] == out$summary["alpha",3])) +# test jags.View--------------------------------------------------------------- +at_home <- identical( Sys.getenv("AT_HOME"), "TRUE" ) +if(at_home){ + ref <- readRDS("longley_reference_fit.Rds") + test <- jags.View(ref) + expect_equal(ncol(test), 10) + expect_equal(colnames(test)[7:10], c("overlap0", "f", "Rhat", "n.eff")) + test2 <- jags.View(out) + expect_equal(ncol(test2), 8) + expect_equal(colnames(test)[7:8], c("overlap0", "f")) +} + # Error when user tries to set seed-------------------------------------------- expect_error(jags(data = data, inits = inits, parameters.to.save = params, model.file = modfile, n.chains = 1, n.adapt = 100, n.iter = 100, @@ -244,8 +256,6 @@ expect_equal(nrow(out$summary), 21) expect_equal(rownames(out$summary)[20], "mu[1,2]") # When no stochastic nodes (and deviance is not calculated)-------------------- -library(jagsUI) - data <- list(a = 1, b = Inf) modfile <- tempfile() diff --git a/inst/tinytest/test_plots.R b/inst/tinytest/test_plots.R new file mode 100644 index 0000000..4e06f39 --- /dev/null +++ b/inst/tinytest/test_plots.R @@ -0,0 +1,49 @@ +# Handle missing values in plots----------------------------------------------- +ref <- readRDS('longley_reference_fit.Rds') + +# Rhat_min +pdf(NULL) +test1 <- traceplot(ref, Rhat_min=1.0005) +expect_equal(test1, NULL) +expect_error(traceplot(ref, Rhat_min=1.1)) +dev.off() + +# Plotting array parameter +pdf(NULL) +test2 <- traceplot(ref, "mu", ask=FALSE) +expect_equal(test2, NULL) +dev.off() + +# Plotting range +pdf(NULL) +test3 <- traceplot(ref, "mu[1:4]", ask=FALSE) +expect_equal(test3, NULL) +dev.off() + +# Handling NAs +for (i in 1:3){ + ref$samples[[i]][,'alpha'] <- NA +} +ref$samples[[1]][,"beta"] <- NA +ref$samples[[1]][1,"sigma"] <- NA +ref$summary["alpha","Rhat"] <- NA + +pdf(NULL) +test4 <- traceplot(ref, c("alpha", "beta", "sigma")) +expect_equal(test4, NULL) +test5 <- whiskerplot(ref, c("alpha", "beta", "sigma")) +expect_equal(test5, NULL) +# Not ideal, maybe fix someday +test6 <- densityplot(ref, c("alpha", "beta", "sigma")) +expect_equal(test6, NULL) +test7 <- plot(ref, ask=FALSE) +expect_equal(test7, NULL) +dev.off() + +# Incorrect object +expect_error(traceplot(ref$samples)) +# Parameter not found +expect_error(traceplot(ref, "fake")) +expect_error(densityplot(ref, "fake")) +expect_error(whiskerplot(ref, "fake")) +expect_error(plot(ref, "fake")) |