rstan で回帰 連続変数のカテゴリ化バージョン

rstan で回帰　連続変数のカテゴリ化バージョン

data {
  int N;
  int K;
  int C;
  // int G;
  matrix[N,K] X;
  int y[N];
  int start_flg[K];
  int order_flg[K];
  
  //int group[N];
  int Level_no;
  //number of room
  int m[N];
  
}

parameters{
  
  real a;
  vector<lower=-5,upper=5>[K-C] b;
  //real<lower=0> sigma;
  real<lower=0,upper=0.2> sigma_f;
  //real<lower=0,upper=0.5> lambda_f;
  real<lower=0> sigma3;
  vector<lower=0,upper=5>[C] sigma_c;
  vector<lower=0,upper=30>[N] rrr;
  real<lower=0.0,upper=1.0> q;
  //real<lower=0,upper=1>[N] p0;
  vector[N] f;
  //vector[G] f_c;
}

transformed parameters{
  
  vector[N] mu;
  vector[K] c;
  vector<lower=0>[K-C] sigma_b;
  real p0;
  real<lower=0,upper=1> pp0b;
  vector<lower=0,upper=1>[N] pp0;
  //matrix<lower=0,upper=1>[Level_no,N] theta;
  simplex[Level_no] theta[N];
  
  for(i in 1:K){
    if(order_flg[i]==0){
      
      if(start_flg[i]!=0){
        c[i]=0;
      }else{
        c[i]=b[i-sum(start_flg[1:i])];
        sigma_b[i-sum(start_flg[1:i])]=sigma_c[sum(start_flg[1:i])];
      }
    }else{
      
      if(start_flg[i]!=0){
        c[i]=0;
      }else{
        c[i]=c[i-1]+b[i-sum(start_flg[1:i])];
        sigma_b[i-sum(start_flg[1:i])]=sigma_c[sum(start_flg[1:i])];
      }
      
    }
  }
  //c[K]=-sum(b[1:(K-1)]);
  //for(i in 1:N){
    //  f[i]=f_c[group[i]];
    //
      //}
  
  
  mu=X*c+a+1*f;
  
  
  for(i in 1:N){
    //print(i);
    //print(f[i]);
    //print(mu[i]);
    //pp0[i]=inv_logit(mu[i]);
    pp0[i]=1.0/(1+exp(-mu[i]));
    pp0b=1-(1-pp0[i])^rrr[i];
    //pp0b=1.0/(1+exp(-mu[i]));
    
    //print("pp0b");
    //print(pp0b);
    p0=binomial_lpmf(1|m[i],pp0b);
    // print("exponential");
    //print(exp(p0));
    //print("rr");
    //print(rrr[i]);
    //print("q");
    //print(q);
    //theta[i][1]=q+(1-q)*(1-(1-exp(p0))^(rrr[i]));
    theta[i][1]=q+(1-q)*(exp(p0)/(1-exp(binomial_lpmf(0|m[i],pp0b))));
    //print("theta");
    // print(theta[i][j]);
    //pp0[i]=inv_logit(mu[i]);
    
  }
  
  
  for(i in 1:N){
    for(j in 2:Level_no){
      if(j<=m[i]){
        //pp0[i]=inv_logit(mu[i]);
        pp0[i]=1.0/(1+exp(-mu[i]));
        //pp0b=pp0[i];
        pp0b=1-(1-pp0[i])^rrr[i];
        //  print("pp0b");
        //print(pp0b);
        p0=binomial_lpmf(j|m[i],pp0b);
        //print("ij");
        //print(i);
        //print(j);
        //print("mm");
        //print(m[i]);
        //print("rr")
        //print(rrr[i]);
        //print("p0");
        //print(exp(p0));
        //if(p0>=-9){
          theta[i][j]=(1-q)*(exp(p0)/(1-exp(binomial_lpmf(0|m[i],pp0b))));
          //}else{
            //   print("PPP000");
            //   theta[i][j]=(1-q)*0.0001;
            
            //}
      }else{
        theta[i][j]=0;
        
      }
      //print("theta");
      // print(theta[j,i]);
      
    }
  }
  
}

model{
  sigma3 ~ student_t(4,0,3);
  sigma_c ~ normal(0,sigma3);
  rrr ~ gamma(5.6,1.0);
  //rrr ~ normal(4,0.00001);
  //sigma_f ~ exponential(lambda_f);
  b ~ double_exponential(0,sigma_b);
  //b ~ double_exponential(0,1);
  //b ~ double_exponential(0,0.3);
  //b ~ student_t(4,0,0.5);
  //sigma_f ~ student_t(4,0,0.3);
  sigma_f ~ student_t(4,0,3);
  
  //for(i in 1:(K-1)){
    //b ~ normal(0,sigma2);
    //}
  //b ~ rnorm(0,sigma2);
  //f_c ~ normal(0,sigma_f);
  f ~ normal(0,sigma_f);
  //y ~ neg_binomial_2_log(mu, 1.0);
  for(i in 1:N){
    //print("smaple");
    //print(theta[i]);
    y[i] ~ categorical(theta[i]);
    //print(y[i]);
  }
}