data{ // num_obs_total: number of trials int num_obs_total ; // obs: observation on each trial vector[num_obs_total] obs ; // num_subj: number of subj int num_subj ; // num_rows_uW: num rows in uW int num_rows_uW ; // num_cols_uW: num cols in uW int num_cols_uW ; // uW: unique entries in the within predictor matrix matrix[num_rows_uW,num_cols_uW] uW ; // index: index of each trial in flattened subject-by-condition value matrix int obs_index[num_obs_total] ; } transformed data{ // obs_mean: mean obs value real obs_mean = mean(obs) ; // obs_sd: sd of obss real obs_sd = sd(obs) ; // obs_: observations scaled to have zero mean and unit variance vector[num_obs_total] obs_ = (obs-obs_mean)/obs_sd ; //uW_rows_repeated_per_subj: an array where each element is a row from uW repeated for num_subj rows matrix[num_subj,num_cols_uW] uW_rows_repeated_per_subj[num_rows_uW] ; for(row in 1:num_rows_uW){ uW_rows_repeated_per_subj[row] = rep_matrix(uW[row],num_subj) ; } } parameters{ //for parameters below, trailing underscore denotes that they need to be un-scaled in generated quantities // mean_coef_: mean (across subj) for each coefficient row_vector[num_cols_uW] mean_coef_ ; // sd_coef_: sd (across subj) for each coefficient vector[num_cols_uW] sd_coef_ ; // chol_corr: population-level correlations (on cholesky factor scale) amongst within-subject predictors cholesky_factor_corr[num_cols_uW] chol_corr ; // multi_normal_helper: a helper variable for implementing non-centered parameterization matrix[num_cols_uW,num_subj] multi_normal_helper ; // noise_: measurement noise real noise_ ; } transformed parameters{ // compute coefficients for each subject/condition matrix[num_subj,num_cols_uW] subj_coef_ = ( rep_matrix(mean_coef_,num_subj) + transpose( diag_pre_multiply(sd_coef_,chol_corr) * multi_normal_helper ) ) ; matrix[num_rows_uW,num_subj] value_for_subj_cond1 ; profile("v1"){ // Loop over subj and conditions to compute unique entries in design matrix for(this_subj in 1:num_subj){ for(this_condition in 1:num_rows_uW){ value_for_subj_cond1[this_condition,this_subj] = dot_product( subj_coef_[this_subj] , uW[this_condition] ) ; } } } matrix[num_rows_uW,num_subj] value_for_subj_cond2 ; profile("v2"){ // Loop over subjs to compute unique entries in design matrix for(this_subj in 1:num_subj){ value_for_subj_cond2[,this_subj] = rows_dot_product( rep_matrix( subj_coef_[this_subj] , num_rows_uW ) , uW ) ; } } matrix[num_rows_uW,num_subj] value_for_subj_cond3 ; profile("v3"){ // Loop over conditions to compute unique entries in design matrix for(this_condition in 1:num_rows_uW){ value_for_subj_cond3[this_condition] = transpose( rows_dot_product( subj_coef_ , uW_rows_repeated_per_subj[this_condition] ) ); } } } model{ //// // Priors //// // low-near-zero prior on measurement noise noise_ ~ weibull(2,1) ; // weibull(2,1) is peaked around .8 // independent (unpooled) normal(0,1) priors on all elements of mean_coef mean_coef_ ~ std_normal() ; // independent (unpooled) normal(0,1) priors on all elements of sd_coef sd_coef_ ~ std_normal() ; // relatively flat prior on correlations chol_corr ~ lkj_corr_cholesky(2) ; // multi_normal_helper must have normal(0,1) prior for non-centered parameterization to_vector(multi_normal_helper) ~ std_normal() ; // Likelihood obs_ ~ normal( to_vector(value_for_subj_cond1)[obs_index] , noise_ ) ; } generated quantities{ // cor: correlation matrix for the full set of within-subject predictors corr_matrix[num_cols_uW] cor = multiply_lower_tri_self_transpose(chol_corr) ; // sd_coef_: sd (across subj) for each coefficient vector[num_cols_uW] sd_coef = sd_coef_ * obs_sd ; // coef_mean: mean (across subj) for each coefficient row_vector[num_cols_uW] mean_coef = mean_coef_ * obs_sd ; mean_coef[1] = mean_coef[1] + obs_mean ; //adding the intercept (assumes contrast matrix had an intercept column) // noise: measurement noise real noise = noise_ * obs_sd ; // tweak cor to avoid rhat false-alarm for(i in 1:num_cols_uW){ cor[i,i] += uniform_rng(1e-16, 1e-15) ; } }