make vignette for different types of prior value generator (can be inserted into workflow.rmd in the future)

[hyunjimoon]

Using SBC library structure, write workflow where prior value simulation is performed given hyperparameter and only its value such as below. Refer to this working code to adjust the input format.

customprival_generator <- function(true_hp, n_datasets = n_datasets){
  dist2samp <- function(true_hp){
    estBetaParams <- function(mu, var) {
      alpha <- ((1 - mu) / var - 1 / mu) * mu ^ 2
      beta <- alpha * (1 / mu - 1)
      return(params = list(alpha = alpha, beta = beta))
    }
    p21 = rnorm(1, true_mean_sd["mean"][1])  #rnorm(1, true_hp[[1]][1], true_hp[[1]][2]) #rbeta(1,p21_hp$alpha, p21_hp$beta)
    p31 = rnorm(1, true_hp[[2]][1], true_hp[[2]][2])  #rbeta(1, p31_hp$alpha, p31_hp$beta)
    for (i in 1:3){
      rate[i] = rnorm(1,true_hp[[3]][1], true_hp[[3]][2]/3)  #rgamma(1, ratet_hp$shape, ratet_hp$scale)
    }
    return(list(p21, p31, rate))
  }
  #lapply(1:5, function(x) {dist2samp(true_hp)
  true_samples <- list()
  for (i in 1:n_datasets){
    true_samples <- c(true_samples, list(dist2samp(true_hp)))
  }
  true_samples
}
true_samples <- customprival_generator(true_mean_sd, 10)

hp_generator <- function(iter, true_samples){
    # Metadata
    N = 3069
    T = 31
    S = 3
    P = 4
    obs2time = complete(generateMice(), 1)$age_ind
    initial_state = 1
    init_state_vec = rep(0, S);
    init_state_vec[initial_state] = 1
    latent_states = array(NA, dim=c(T,S))
    for (i in n_datasets){
      p21 = true_samples[[i]][[1]] #rbeta(1,p21_hp$alpha, p21_hp$beta)
      p31 = rnorm(1, true_hp[[2]][1], true_hp[[2]][2] )  #rbeta(1, p31_hp$alpha, p31_hp$beta)
      for (i in 1:3){
        #ratet_hp <- gammaParamsConvert(mean=true_hp[[3]][1],sd=true_hp[[3]][2])
        rate[i] = true_hp[[3]][i]  #rnorm(1,true_hp[[3]][1], true_hp[[3]][2]/3)  #rgamma(1, ratet_hp$shape, ratet_hp$scale)
      }
      # Maintenance
      Mnt = array(c(1, 0, 0, p21, 1- p21, 0, p31, 1-p31,0), dim=c(S,S))
      # Deterioration
      tmp_p <- array(rep(0, S * S), dim=c(S,S))
      tmp_p[1,1] = exp(-rate[1]- rate[2])
      tmp_p[2,1] = rate[1] * exp(-rate[3]) * (1-exp(-(rate[1]+ rate[2] - rate[3]))) / (rate[1]+ rate[2] - rate[3])
      tmp_p[3,1] = exp(-rate[3])
      Dtr <- array(c(tmp_p[1,1], tmp_p[2,1], 1- tmp_p[1,1], 0, tmp_p[3,1], 1 - tmp_p[3,1], 0,0,1), dim=c(S, S))
      latent_states[1,] = Dtr %*% init_state_vec
      for (t in 2:T){
        latent_states[t,] =  (Dtr %*% Mnt) %*% latent_states[t-1,]
      }
    }
    #p21_hp <- estBetaParams(true_hp[[1]][1],(true_hp[[1]][2])^2) this gives all 1
    #p31_hp <- estBetaParams(true_hp[[2]][1],(true_hp[[2]][2])^2)

    list(
      generated = lapply(seq(1:length(obs2time)), function(n){sample(c(1,2,3), 1, prob = latent_states[obs2time[n],])}),
      parameters = list(
        p21 = p21,
        p31 = p31,
        rate = rate
      )
    )
  }