Skip to content
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 20 additions & 0 deletions R/calib.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
library(cmdstanr)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

instead of importing libraries, try to use namespace::func_name.
Additionally, include necessary libraries into Imports in DESCRIPTION and add ROxygen documentation.

library(dplyr)
library(tidyverse)
library(reshape)

iter_gen_inf <- function(sbc_obj, prior_theta, par_names, N, M, data){ #fixed params as prior_phi and add clampedStan
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we can just receive sampled_y and post_theta as arguments instead of SBCObject and re-fitting.

sampled_y = sbc_obj$sample_y_tilde(prior_theta, data=data)
post_theta = sbc_obj$sample_theta_bar_y(sampled_y, data=data, pars=par_names, fit_iter=M)
NMP_tb <- as_tibble(post_theta)
NMP_G <- NMP_tb %>% gather(variable)
NMP_G$par <- gsub("\\..*","",NMP_G$variable)
NMP_G$prior <- sapply(NMP_G$variable, function(x){prior_theta[as.integer(gsub(".*\\.","",x)), gsub("\\..*","",x)]})
Epost_bar_pri <- NMP_G %>% group_by(par, prior) %>% summarise('mean' = round(mean(value),3))
Varpost_bar_pri<- NMP_G %>% group_by(par, prior) %>% summarise('sd' = sd(value))
postVar <- inner_join(Epost_bar_pri %>% group_by(par) %>% summarise('VE' = sd(mean)), Varpost_bar_pri %>% group_by(par) %>% summarise('EV' = mean(sd^2))) %>%
mutate("postVar" = VE + EV) %>% select(par, postVar)
priVar = NMP_G %>% select(prior, par) %>% group_by(par) %>% summarise('priVar' = sd(prior)^2)
ratioVar = inner_join(postVar, priVar) %>% mutate("ratioVar" = postVar/priVar)
return(ratioVar)
}