# SBC Interpretation:bias of prior and bias of MCMC sampling are trade off

**URL:** <https://discourse.mc-stan.org/t/sbc-interpretation-bias-of-prior-and-bias-of-mcmc-sampling-are-trade-off/9712>\
**Category:** Modeling\
**Tags:** simulation-based-calibration\
**Created:** [July 14, 2019, 9:07am UTC](https://discourse.mc-stan.org/t/sbc-interpretation-bias-of-prior-and-bias-of-mcmc-sampling-are-trade-off/9712 "2019-07-14T09:07:28Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jean\_Billie](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jean_billie/32/9186_2.png) [@Jean\_Billie](https://discourse.mc-stan.org/u/Jean_Billie)\
**Post date:** [July 14, 2019, 9:07am UTC](https://discourse.mc-stan.org/t/sbc-interpretation-bias-of-prior-and-bias-of-mcmc-sampling-are-trade-off/9712/1 "2019-07-14T09:07:28Z")

</div>

I implemented SBC for my model.  
The following is the result. Please tell me is there any problem? In the MCMC simulation, the divergent transition occurred.

In my model, when I fit my model to data, then the strong assumption on variance of prior is not required. However, when I implement SBC, if I put large variances on prior, then it generated a very odd data, which cause the very odd sampling. So, in my model, to implement SBC, I have to carefully choose the prior. My opinion, SBC algorithm requires the suitable informative priors and SBC is strongly affected the selection of prior.

If I choose informative prior, then SBC tells me it is good MCMC sampling, so, sampling itself contains no bias, however, such informative prior itself contains bias. And if I do not assume informative prior, then the odd data are made by the noninoformative priors and MCMC sampling contains bias. So, I think bias of prior and bias of MCMC sampling is trade off.

 ![Rplot01](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/c/c9805ad08232301a116a0e08c419238df8bdfce1.jpeg) R code

```
stanmodel <- stan_model("sbc.stan")

tttttt <- function( ww=-0.81,www =0.001, #0.001 can be change to 0.1 but 1 is too large to get sampling.
                 mm=0.65,mmm=0.001,
                 vv=5.31,vvv=0.001,
                 zz= 1.55,zzz=0.001 ){
  
  
  output <- sbc(stanmodel, data = list(
    www= www,
    mmm= mmm,
    vvv= vvv,
    zzz =zzz,
    
    ww=ww,
    mm=mm,
    vv=vv,
    zz=zz,
    
    N = 3, NL = 259, NI = 57,C=3,c=c(3,2,1)), M = 500, refresh = 0)
}

fit <- tttttt( ww=-0.81,www =0.001,
                 mm=0.65,mmm=0.001,
                 vv=5.31,vvv=0.001,
                 zz= 1.55,zzz=0.001 )

plot(fit, bins = 111)

```

* * *

**sbc.stan**

```
   data{ // SBC

    //This is not prior truth data, but somedata to run
      int <lower=0>N; //This is exactly same as C
      int <lower=0>NL; //Number of Binomial trials
      int <lower=0>NI; //This is redandunt 
      int <lower=0>C; // Number of Confidence level
      int <lower=0>c[N]; //Each component means confidence level

    //Prior which shold be specified
    real www;
    real mmm;
    real vvv;
    real zzz;
    real zz;
    real ww;
    real vv;
    real mm;
     
    }

    transformed data {
     int h[C];
     int f[C];

      real w_ ; 
      real <lower=0>dz_[C-1] ; 
      real m_; 
      real <lower =0> v_;  

      real <lower=0,upper=1>p_[C];
      real <lower=0>l_[C];
      real <lower=0>dl_[C];
       real z_[C];
     

    real a_;
    real <lower=0>b_;

                              w_ = normal_rng (ww,www);
        for(cd in 1:C-1) dz_[cd] = normal_rng (zz,zzz);
                              m_ = normal_rng (mm,mmm);
                              v_ = normal_rng (vv,vvv);

    a_=m_/v_;
    b_=1/v_;

         for(cd in 1 : C-1) { z_[1]=w_;
                          z_[cd+1] =z_[cd] +dz_[cd];
                          }

         for(cd in 1 : C) { if (cd==C) {
                               p_[cd] = 1 - Phi((z_[cd] - m_)/v_);
                                          }else{
                               p_[cd] = Phi((z_[cd+1] - m_)/v_)- Phi( (z_[cd] -m_)/v_);

                                               }
                           }

         for(cd in 1 : C) {l_[cd] = (-1)*log(Phi(z_[cd])); }
         for(cd in 1:C){
                     if (cd==C) {dl_[cd]=fabs(l_[cd]-0);
                     }else{

           dl_[cd]=fabs(l_[cd]-l_[cd+1]);
              }
              }

               for(n in 1:N) {
        h[n] = binomial_rng(NL, p_[c[n]]);
     // fff[n] ~ poisson( l[c[n]]*NL);//Non-Chakraborty's model
        f[n] = poisson_rng (dl_[c[n]]*NI);//Chakraborty's model //<-------very very very coution, not n but c[n] 2019 Jun 21
     // fff[n] ~ poisson( l[c[n]]*NI);//Non-Chakraborty's model

                           }
    }

    parameters{
      real w;
      real <lower =0>dz[C-1];
      real m;  
      real <lower=0>v;

    }

    transformed parameters {

      real <lower=0,upper=1>p[C];
      real <lower=0>l[C];
      real <lower=0>dl[C];
      real z[C];

    real a;
    real b;

    a=m/v;
    b=1/v;

         for(cd in 1 : C-1) { z[1] = w;
                             z[cd+1] = z[cd] +dz[cd];
                            }

         for(cd in 1 : C) {
           if (cd==C) { p[cd] = 1 - Phi((z[cd] -m)/v);
         }else{
                               p[cd] = Phi((z[cd+1] -m)/v)- Phi((z[cd] -m)/v);

         }
         }

         for(cd in 1 : C) { l[cd] = (-1)*log(Phi(z[cd])); }
         for(cd in 1:C){
                  if (cd==C) {dl[cd] = fabs(l[cd]-0);
                     }else{

                              dl[cd] = fabs(l[cd]-l[cd+1]);
              }
              }
    }

    model{
           for(n in 1:N) {
                             h[n] ~ binomial(NL, p[c[n]]);
     // fff[n] ~ poisson( l[c[n]]*NL);//Non-Chakraborty's model
                             f[n] ~ poisson(dl[c[n]]*NI);//Chakraborty's model //<-------very very very coution, not n but c[n] 2019 Jun 21
     // fff[n] ~ poisson( l[c[n]]*NI);//Non-Chakraborty's model

                           }

       // priors  

                              w ~ normal(ww,www);
        for(cd in 1:C-1) dz[cd] ~ normal(zz,zzz);
                              m ~ normal(mm,mmm);
                              v ~ normal(vv,vvv);

    }

    generated quantities { // these adhere to the conventions above
    int h_[C];
    int f_[C];
    vector [1] pars_;
    int ranks_[1] = {w > w_};
    pars_[1] = w_;
    h_ = h;
    f_ = f;
     
    }
```