Skip to content

self-consistent region iteration#15

Open
hyunjimoon wants to merge 1 commit intomasterfrom
selfCalib
Open

self-consistent region iteration#15
hyunjimoon wants to merge 1 commit intomasterfrom
selfCalib

Conversation

@hyunjimoon
Copy link
Owner

@hyunjimoon hyunjimoon commented May 6, 2021

What? find the self-consistency parmater regions via contraction mapping
Assumption: variance of parameter samples would be similar before and after in self-consistency parmater regions
How? run until the ratio < 1.1 for all parameters

Test code:

model_ncp = cmdstanr::cmdstan_model("tests/eightschools_ncp.stan")
sbc_obj_ncp = SBC::SBCModel$new(name="eightschools_ncp", stan_model=model_ncp)

model_cp = cmdstanr::cmdstan_model("tests/eightschools_cp.stan")
sbc_obj_cp = SBC::SBCModel$new(name="eightschools_cp", stan_model=model_cp)

par_names <- list("theta", "mu", "tau") #unlist(dimnames(prior_theta)[[2]]) #pars = list("theta", "mu", "tau")

N = 100
M = 10

prior_theta_ncp = sbc_obj_ncp$sample_theta_tilde_stan(par_names, N, data=data) # partial pars possible
ncp_ratio <- iter_gen_inf(sbc_obj_ncp, prior_theta_ncp, par_names, N, M, data)

prior_theta_cp = sbc_obj_cp$sample_theta_tilde_stan(par_names, N, data=data) # partial pars possible
cp_ratio <- iter_gen_inf(sbc_obj_cp, prior_theta_cp, par_names, N, M, data)
  1. N,M = (10,10)
> ncp_ratio
# A tibble: 10 x 4
   par      postVar priVar ratioVar
   <chr>      <dbl>  <dbl>    <dbl>
 1 mu          9.22   39.5    0.233
 2 tau        38.0   126.     0.302
 3 theta[1]   66.5   111.     0.601
 4 theta[2]   81.4   721.     0.113
 5 theta[3]   76.2   102.     0.744
 6 theta[4]   39.9    37.9    1.05 
 7 theta[5]   42.9    50.0    0.860
 8 theta[6]   87.9    37.9    2.32 
 9 theta[7]   72.7    72.2    1.01 
10 theta[8]  112.    628.     0.178

> cp_ratio
# A tibble: 10 x 4
   par      postVar priVar ratioVar
   <chr>      <dbl>  <dbl>    <dbl>
 1 mu          11.2   22.5   0.498 
 2 tau         58.4  414.    0.141 
 3 theta[1]    56.5  971.    0.0582
 4 theta[2]    89.9 4240.    0.0212
 5 theta[3]    79.8 2056.    0.0388
 6 theta[4]    51.1  221.    0.232 
 7 theta[5]    30.4  156.    0.195 
 8 theta[6]    32.8   97.6   0.336 
 9 theta[7]    33.9  194.    0.174 
10 theta[8]    35.9   84.6   0.424 
  1. N,M = (100,10)
> ncp_ratio
# A tibble: 10 x 4
   par      postVar   priVar ratioVar
   <chr>      <dbl>    <dbl>    <dbl>
 1 mu          8.57     20.6 0.415   
 2 tau       119.    39191.  0.00304 
 3 theta[1]  162.     8041.  0.0201  
 4 theta[2]  119.     6061.  0.0196  
 5 theta[3]  171.     8562.  0.0199  
 6 theta[4]  146.     2253.  0.0650  
 7 theta[5]  137.    83883.  0.00163 
 8 theta[6]  141.   170061.  0.000829
 9 theta[7]   89.4  121873.  0.000734
10 theta[8]  148.     5768.  0.0257  
> cp_ratio
# A tibble: 10 x 4
   par      postVar priVar ratioVar
   <chr>      <dbl>  <dbl>    <dbl>
 1 mu          12.4   22.7   0.544 
 2 tau         62.2 3795.    0.0164
 3 theta[1]    68.7 2550.    0.0270
 4 theta[2]    72.2 1269.    0.0569
 5 theta[3]    65.7 2647.    0.0248
 6 theta[4]   112.  4385.    0.0255
 7 theta[5]    94.3 1618.    0.0583
 8 theta[6]    91.4 4337.    0.0211
 9 theta[7]    58.2  968.    0.0601
10 theta[8]    89.2 4971.    0.0180

@hyunjimoon hyunjimoon requested a review from Dashadower May 6, 2021 05:05
@@ -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(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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants