/* History:
2025-10-30 Add DEBUG, check the existence of files loaded, and add code to call modules in
           source_v143.sas for alternative v14.3 code.
		   Check the log argument in marginal_likelihood.
*/

/*
Program Purpose:
IML modules for the main MCMC loop
*/

proc iml;

	**Main Gibbs sampler loop to fit multivar MCMC model;
	start mcmc_main_loop(outname,
											 beta_mean_prior,
											 beta_covariance_prior,
											 r_matrix_prior,
											 theta_matrix_prior,
											 v_matrix_prior,
											 sigma_u_prior,
											 sigma_e_prior,
											 u_matrix_prior,
											 w_matrix_prior,
											 xbeta_prior,
											 xbeta_u_prior,
											 recall_availability,
											 subject_weighting,
											 episodic_indicator_matrices,
											 covariate_matrices,
											 wt_covariate_matrices,
											 wt_covariate_sq_sums,
											 num_mcmc_iterations,
											 num_burn,
											 num_thin,
											 num_post,
											 num_subjects,
											 num_recalls,
											 num_episodic,
											 num_daily,
											 num_variables,
											 num_coefficients,
											 sigma_u_constant,
											 save_mcmc_u,
											 save_all_u,
											 do_log_likelihood,
											 variable_names,
											 coefficient_names,
											 has_never_consumers,
											 alpha1_mean_prior,
											 alpha1_covariance_prior,
											 consumer_probabilities_prior,
											 never_consumer_covariate_matrix,
											 num_never_consumer_coefficients,
											 never_consumer_coefficient_names,
											 low_memory_mode);
	
		if symgetn("DEBUG") then print "begin module: mcmc_main_loop";			
		alg = symgetn("v143code") + 1;
	
		**Initialize parameters to their priors;
		beta = beta_mean_prior;
		
		if num_episodic > 1 then do;
		
			r_matrix = r_matrix_prior;
			theta_matrix = theta_matrix_prior;
		end;
		else do;
		
			r_matrix = {};
			theta_matrix = {};
		end;
		
		v_matrix = v_matrix_prior;
		sigma_e = sigma_e_prior;
		sigma_u = sigma_u_prior;
		u_matrix = u_matrix_prior;
		w_matrix = w_matrix_prior;
		xbeta = xbeta_prior;
		xbeta_u = xbeta_u_prior;
		
		thinned_iterations = do(num_burn + 1, num_mcmc_iterations, num_thin);
		num_thinned_iterations = int((num_mcmc_iterations - (num_burn + 1))/num_thin) + 1;
		
		**Initialize parameter trace data;
		
		**beta;
		beta_mean = ListCreate(num_variables);
		beta_variance = ListCreate(num_variables);
		
		beta_trace = ListCreate(num_variables);
		beta_trace_names = ListCreate(num_variables);
		
		do var_j = 1 to num_variables;
		
			beta_mean$var_j = j(num_coefficients[var_j], 1, 0);
			beta_variance$var_j = j(num_coefficients[var_j], 1, 0);
			
			beta_j_trace_names = j(num_coefficients[var_j], 1, BlankStr(32));
			
			coefficient_names_j = coefficient_names$var_j;
			
			intercept_j = (coefficient_names_j[1] = "intercept");
			num_covariates_j = num_coefficients[var_j] - intercept_j;
			
			if intercept_j = 1 then beta_j_trace_names[1] = cats("beta", var_j, "_intercept");
			if num_covariates_j > 0 then beta_j_trace_names[intercept_j + (1:num_covariates_j)] = cats("beta", var_j, "_covariate", 1:num_covariates_j);
			
			if upcase(low_memory_mode) = "Y" then do;
			
				beta_trace$var_j = cats(outname, "_beta", var_j);
				create (beta_trace$var_j) var (beta_j_trace_names);
			end;
			else beta_trace$var_j = j(num_mcmc_iterations, num_coefficients[var_j], .);
			
			beta_trace_names$var_j = beta_j_trace_names;
		end;
		
		**sigma-u;
		sigma_u_mean = j(num_variables, num_variables, 0);
		sigma_u_variance = j(num_variables, num_variables, 0);
		
		sigma_u_trace_names = cats("sigma_u_row", symsqr(row(sigma_u)), "_col", symsqr(col(sigma_u)));
		
		if upcase(low_memory_mode) = "Y" then do;
		
			sigma_u_trace = cats(outname, "_sigma_u");
			create (sigma_u_trace) var sigma_u_trace_names;
		end;
		else sigma_u_trace = j(num_mcmc_iterations, num_variables*(num_variables+1)/2);
		
		**sigma-e;
		sigma_e_mean = j(num_variables, num_variables, 0);
		sigma_e_variance = j(num_variables, num_variables, 0);
		
		sigma_e_trace_names = cats("sigma_e_row", symsqr(row(sigma_e)), "_col", symsqr(col(sigma_e)));
		
		if upcase(low_memory_mode) = "Y" then do;
		
			sigma_e_trace = cats(outname, "_sigma_e");
			create (sigma_e_trace) var sigma_e_trace_names;
		end;
		else sigma_e_trace = j(num_mcmc_iterations, num_variables*(num_variables+1)/2);
		
		**Initialize U matrix trace data for main MCMC and post-MCMC iterations;
		if upcase(save_mcmc_u) = "Y" then do;
		
			if upcase(save_all_u) = "Y" then do;
				
				saved_u_traces = 1:num_mcmc_iterations;
				num_mcmc_u = num_mcmc_iterations;
			end;
			else do;
			
				saved_u_traces = thinned_iterations;
				num_mcmc_u = num_thinned_iterations;
			end;
			
			u_mcmc_names = cats("u_col", 1:num_variables);
			
			if upcase(low_memory_mode) = "Y" then do;
			
				u_matrices_mcmc = cats(outname, "_u_mcmc");
				create (u_matrices_mcmc) var ("iteration" || u_mcmc_names);
			end;
			else u_matrices_mcmc = j(num_subjects * num_mcmc_u, num_variables, .);
		end;	
		else u_matrices_mcmc = {};
		
		if num_post > 0 then do;
		
			u_post_names = cats("u_col", 1:num_variables);
		
			if upcase(low_memory_mode) = "Y" then do;
			
				u_matrices_post = cats(outname, "_u_post");
				create (u_matrices_post) var ("post_mcmc_iteration" || u_post_names);
			end;
			else u_matrices_post = j(num_subjects * num_post, num_variables, .);
		end;
		else u_matrices_post = {};
		
		**never-consumer parameters;
		if has_never_consumers = 1 then do;
		
			**set never-consumer parameters to their priors;
			conni1 = {};
			alpha1 = alpha1_mean_prior;
			consumer_probabilities = consumer_probabilities_prior;
			
			**initialize never-consumer parameter trace data;
			
			**alpha1;
			alpha1_mean = j(num_never_consumer_coefficients, 1, 0);
			alpha1_variance = j(num_never_consumer_coefficients, 1, 0);
			
			alpha1_trace_names = j(num_never_consumer_coefficients, 1, BlankStr(32));
			
			never_consumer_intercept = (never_consumer_coefficient_names[1] = "intercept");
			num_never_consumer_covariates = num_never_consumer_coefficients - never_consumer_intercept;
			
			if never_consumer_intercept then alpha1_trace_names = cats("alpha1_intercept");
			if num_never_consumer_covariates > 0 then alpha1_trace_names[never_consumer_intercept + (1:num_never_consumer_covariates)] = cats("alpha1_covariate", 1:num_never_consumer_covariates);
			
			if upcase(low_memory_mode) = "Y" then do;
			
				alpha1_trace = cats(outname, "_alpha1");
				create (alpha1_trace) var alpha1_trace_names;
			end;
			else alpha1_trace = j(num_mcmc_iterations, num_never_consumer_coefficients, .);
			
			**consumer probabilities;
			consumer_probabilities_mean = 0;
			consumer_probabilities_variance = 0;
			
			conprob1_trace_name = "conprob1";
			
			if upcase(low_memory_mode) = "Y" then do;
			
				consumer_probabilities_trace = cats(outname, "_conprob1");
				create (consumer_probabilities_trace) var conprob1_trace_name;
			end;
			else consumer_probabilities_trace = j(num_mcmc_iterations, 1, .);
		end;
		else do;
		
			conni1 = {};

			alpha1 = {};
			alpha1_mean = {};
			alpha1_variance = {};
			alpha1_trace = {};
			
			consumer_probabilities = {};
			consumer_probabilities_mean = {};
			consuemr_probabilities_variance = {};
			consumer_probabilities_trace = {};
		end;
		
		**Gibbs sampler loop;	
		thinned_iter = 0;
			 
		do iter = 1 to num_mcmc_iterations + num_post;
		
			if has_never_consumers = 1 then do;
			
							**Update Ni;
							conni1 = update_conni1(alpha1,
														 				 consumer_probabilities,
														 				 xbeta_u,
														 				 never_consumer_covariate_matrix,
														 				 episodic_indicator_matrices,
														 				 recall_availability,
														 				 num_subjects,
														 				 num_recalls);
														 				 
							**Update alpha1;
							alpha1 = update_alpha1(alpha1_mean_prior,
														 				 alpha1_covariance_prior,
														 				 conni1,
														 				 subject_weighting,
														 				 never_consumer_covariate_matrix,
														 				 num_never_consumer_coefficients);
														 				 
							**Update consumer probabilities;
							consumer_probabilities = update_consumer_probabilities(never_consumer_covariate_matrix,
																														 				 alpha1);
			end;

			**Update W matrix;
			w_matrix = update_w_matrix(w_matrix,
																 xbeta_u,
																 sigma_e,
																 recall_availability,
																 episodic_indicator_matrices,
																 num_subjects,
																 num_episodic,
																 num_variables,
																 num_recalls);
			
			if iter <= num_mcmc_iterations then do;
			
				**Calculate W-XBeta-U cross-residual sum;
				**The W-XBeta-U calculation corresponds to the error (epsilon) term in Equation 3.5 of Zhang, et al. (2011);

				w_cross_residual_sum = calculate_w_cross_residual_sum(w_matrix,
																															xbeta_u,
																															recall_availability,
																															subject_weighting,
																															num_variables,
																															num_recalls);
				
				**Update r matrix;
				r_matrix = update_r_matrix(r_matrix,
																	 theta_matrix,
																	 v_matrix,
																	 w_cross_residual_sum,
																	 recall_availability,
																	 subject_weighting);
				
				**Update theta matrix;
				theta_matrix = update_theta_matrix(theta_matrix,
																					 r_matrix,
																					 v_matrix,
																					 w_cross_residual_sum);
				
				**Update v matrix;
				v_matrix = update_v_matrix(v_matrix,
																	 r_matrix,
																	 theta_matrix,
																	 w_cross_residual_sum,
																	 recall_availability,
																	 subject_weighting,
																	 num_episodic,
																	 num_daily);
				
				**Update sigma-e;
/* BW CHANGE 2025-10-16 */
				*sigma_e = v_matrix * v_matrix`;
				sigma_e = matMultAAt(v_matrix, alg);
				
				**Update sigma-u;
				sigma_u = update_sigma_u(sigma_u,
																 sigma_u_prior,
																 u_matrix,
																 subject_weighting,
																 num_subjects,
																 num_variables,
																 sigma_u_constant);

			end;
			
			**Update U matrix;
			u_matrix = update_u_matrix(u_matrix,
																 sigma_u,
																 sigma_e,
																 w_matrix,
																 xbeta,
																 recall_availability,
																 num_subjects,
																 num_variables,
																 num_recalls,
																 has_never_consumers,
																 conni1);
			
			if iter <= num_mcmc_iterations then do;
			
				**Update beta;
				beta = update_beta(w_matrix,
													 u_matrix,
													 sigma_e,
													 xbeta,
													 beta_mean_prior,
													 beta_covariance_prior,
													 covariate_matrices,
													 wt_covariate_matrices,
													 wt_covariate_sq_sums,
													 recall_availability,
													 num_subjects,
													 num_variables,
													 num_coefficients,
													 num_recalls,
													 has_never_consumers,
													 conni1,
													 beta$1);
			end;
			
			if iter <= num_mcmc_iterations then do;
			
				is_thinned = (iter > num_burn) & (mod(iter, num_thin) = mod(num_burn + 1, num_thin));
				
				if is_thinned = 1 then thinned_iter = thinned_iter + 1;
				
				**Save parameter traces;
				
				**beta;
				do var_j = 1 to num_variables;
				
					beta_j_out = beta$var_j`;
					
					if upcase(low_memory_mode) = "Y" then do;
					
						setout (beta_trace$var_j);
						append from beta_j_out;
					end;
					else do;
					
						beta_j_trace = beta_trace$var_j;
						beta_j_trace[iter,] = beta_j_out;
						beta_trace$var_j = beta_j_trace;
					end;
				end;
				
				**sigma-u;
				sigma_u_out = symsqr(sigma_u)`;
				
				if upcase(low_memory_mode) = "Y" then do;
				
					setout (sigma_u_trace);
					append from sigma_u_out;
				end;
				else sigma_u_trace[iter,] = sigma_u_out;
				
				**sigma-e;
				sigma_e_out = symsqr(sigma_e)`;
				
				if upcase(low_memory_mode) = "Y" then do;
				
					setout (sigma_e_trace);
					append from sigma_e_out;
				end;
				else sigma_e_trace[iter,] = sigma_e_out;
				
				**if specified, save U matrix;
				if upcase(save_mcmc_u) = "Y" then do;
			
					if upcase(save_all_u) = "Y" | is_thinned = 1 then do;
					
						if upcase(low_memory_mode) = "Y" then do;
						
							u_mcmc_out = repeat(iter, num_subjects, 1) || u_matrix;
						
							setout (u_matrices_mcmc);
							append from u_mcmc_out;
							
							free u_mcmc_out;
						end;
						else do;
						
							if upcase("save_all_u") = "Y" then mcmc_u_index = iter;
							else mcmc_u_index = thinned_iter;
							
							u_matrices_mcmc[num_subjects*(mcmc_u_index-1) + (1:num_subjects),] = u_matrix;
						end;
					end;
				end;
				
				**never-consumer parameter traces;
				if has_never_consumers = 1 then do;
				
					**alpha1;
					alpha1_out = alpha1`;
					
					if upcase(low_memory_mode) = "Y" then do;
					
						setout (alpha1_trace);
						append from alpha1_out;
					end;
					else alpha1_trace[iter,] = alpha1_out;
				
					**consumer probabilities;
					consumer_probabilities_all = mean(consumer_probabilities);
				
					if upcase(low_memory_mode) = "Y" then do;
					
						setout (consumer_probabilities_trace);
						append from consumer_probabilities_iter;
					end;
					else consumer_probabilities_trace[iter] = consumer_probabilities_iter;
				end;
				
				**update means and variances at thinned iterations;
				if is_thinned = 1 then do;
					
					**beta;
					do var_j = 1 to num_variables;
					
						beta_j_delta1 = beta$var_j - beta_mean$var_j;
						beta_mean$var_j = beta_mean$var_j + beta_j_delta1/thinned_iter;
						
						beta_j_delta2 = beta$var_j - beta_mean$var_j;
						beta_variance$var_j = beta_variance$var_j + beta_j_delta1 # beta_j_delta2;
					end;
					
					**sigma-u;
					sigma_u_delta1 = sigma_u - sigma_u_mean;
					sigma_u_mean = sigma_u_mean + sigma_u_delta1/thinned_iter;
					
					sigma_u_delta2 = sigma_u - sigma_u_mean;
					sigma_u_variance = sigma_u_variance + sigma_u_delta1 # sigma_u_delta2;
					
					**sigma-e;
					sigma_e_delta1 = sigma_e - sigma_e_mean;
					sigma_e_mean = sigma_e_mean + sigma_e_delta1/thinned_iter;
					
					sigma_e_delta2 = sigma_e - sigma_e_mean;
					sigma_e_variance = sigma_e_variance + sigma_e_delta1 # sigma_e_delta2;
					
					**never-consumer parameters;
					if has_never_consumers = 1 then do;
					
						**alpha1;
						alpha1_delta1 = alpha1 - alpha1_mean;
						alpha1_mean = alpha1_mean + alpha1_delta1/thinned_iter;
						
						alpha1_delta2 = alpha1 - alpha1_mean;
						alpha1_variance = alpha1_variance + alpha1_delta1 # alpha1_delta2;
						
						**consumer probabilities;
						consumer_probabilities_delta1 = consumer_probabilities_all - consumer_probabilities_mean;
						consumer_probabilities_mean = consumer_probabilities_mean + consumer_probabilities_delta1/thinned_iter;
						
						consumer_probabilities_delta2 = consumer_probabilities_all - consumer_probabilities_mean;
						consumer_probabilities_variance = consumer_probabilities_delta1 # consumer_probabilities_delta2;
					end;
				end;
			end;
			else do;
			
				**save U matrix (post-MCMC);
				post_mcmc_iter = iter - num_mcmc_iterations;
				
				if upcase(low_memory_mode) = "Y" then do;
					
					u_post_out = repeat(post_mcmc_iter, num_subjects, 1) || u_matrix;
					
					setout (u_matrices_post);
					append from u_post_out;
					
					free u_post_out;
				end;
				else u_matrices_post[num_subjects*(post_mcmc_iter-1) + (1:num_subjects),] = u_matrix;
			end;
			
			if iter = num_mcmc_iterations then do;
				
				**fix beta, sigma-e, and sigma-u at posterior means;
				beta = beta_mean;
				
				sigma_u = sigma_u_mean;
				
				sigma_e = sigma_e_mean;
				
				**calculate variances;
				do var_j = 1 to num_variables;
				
					beta_variance$var_j = beta_variance$var_j/(num_thinned_iterations - 1);
				end;
				
				sigma_u_variance = sigma_u_variance/(num_thinned_iterations - 1);
				
				sigma_e_variance = sigma_e_variance/(num_thinned_iterations - 1);
				
				if has_never_consumers = 1 then do;
				
					alpha1_variance = alpha1_variance/(num_thinned_iterations - 1);
					
					consumer_probabilities_variance = consumer_probabilities_variance/(num_thinned_iterations - 1);
				end;
			end;
			
			**Update Xbeta;
			if iter <= num_mcmc_iterations then do;
			
				do day_k = 1 to num_recalls;
				
					xbeta_k = xbeta$day_k;
					do var_j = 1 to num_variables;
/* BW CHANGE 2025-10-22 */						
						*xbeta_k[,var_j] = covariate_matrices$var_j$day_k * beta$var_j;
						xbeta_k[,var_j] = matMultAB(covariate_matrices$var_j$day_k, beta$var_j, alg);
					end;
					xbeta$day_k = xbeta_k;
				end;
			end;
			
			**Update Xbeta-U;
			do day_k = 1 to num_recalls;
			
				xbeta_u$day_k = xbeta$day_k + u_matrix;
			end;
		end;
		
		
		**Calculate log-likelihood;
		if upcase(do_log_likelihood) = "Y" then do;
		
			log_likelihood = marginal_likelihood(w_matrix,
																				 	 episodic_indicator_matrices,
																				 	 xbeta,
																				 	 sigma_u,
																				 	 sigma_e,
																				 	 recall_availability,
																				 	 subject_weighting,
																				 	 num_subjects,
																				 	 num_recalls,
																				 	 num_episodic,
																				 	 num_daily,
																				 	 num_variables,
																				 	 has_never_consumers,
																				 	 alpha1_mean,
																				 	 never_consumer_covariate_matrix);
																				 	 
			create (cats(outname, "_ll")) var {log_likelihood};
				append;
			close (cats(outname, "_ll"));
		end;
																			 
		**Output parameters;
		
		**beta;
		do var_j = 1 to num_variables;
		
			_VARIABLE_ = variable_names[var_j];
			coefficient_names_j = coefficient_names$var_j;
		
			beta_j_mean_t = beta_mean$var_j`;
			
			create (cats(outname, "_beta", var_j, "_mean")) from beta_j_mean_t[rowname=_VARIABLE_ colname=coefficient_names_j];
				append from beta_j_mean_t[rowname=_VARIABLE_];
			close (cats(outname, "_beta", var_j, "_mean"));
			
			beta_j_variance_t = beta_variance$var_j`;
			
			create (cats(outname, "_beta", var_j, "_var")) from beta_j_variance_t[rowname=_VARIABLE_ colname=coefficient_names_j];
				append from beta_j_variance_t[rowname=_VARIABLE_];
			close (cats(outname, "_beta", var_j, "_var"));
		
			if upcase(low_memory_mode) = "Y" then close (beta_trace$var_j);
			else do;
			
				beta_j_trace = beta_trace$var_j;
				beta_j_trace_names = beta_trace_names$var_j;
				
				create (cats(outname, "_beta", var_j)) from beta_j_trace[colname=beta_j_trace_names];
					append from beta_j_trace;
				close (cats(outname, "_beta", var_j));
			end;
		end;
		
		**sigma-u;
		_VARIABLE_ = variable_names;
		
		create (cats(outname, "_sigma_u_mean")) from sigma_u_mean[rowname=_VARIABLE_ colname=_VARIABLE_];
			append from sigma_u_mean[rowname=_VARIABLE_];
		close (cats(outname, "_sigma_u_mean"));
		
		create (cats(outname, "_sigma_u_var")) from sigma_u_variance[rowname=_VARIABLE_ colname=_VARIABLE_];
			append from sigma_u_variance[rowname=_VARIABLE_];
		close (cats(outname, "_sigma_u_var"));
		
		if upcase(low_memory_mode) = "Y" then close (sigma_u_trace);
		else do;
		
			create (cats(outname, "_sigma_u")) from sigma_u_trace[colname=sigma_u_trace_names];
				append from sigma_u_trace;
			close (cats(outname, "_sigma_u"));
		end;
		
		**sigma-e;
		_VARIABLE = variable_names;
		
		create (cats(outname, "_sigma_e_mean")) from sigma_e_mean[rowname=_VARIABLE_ colname=_VARIABLE_];
			append from sigma_e_mean[rowname=_VARIABLE_];
		close (cats(outname, "_sigma_e_mean"));
		
		create (cats(outname, "_sigma_e_var")) from sigma_e_variance[rowname=_VARIABLE_ colname=_VARIABLE_];
			append from sigma_e_variance[rowname=_VARIABLE_];
		close (cats(outname, "_sigma_e_var"));
		
		if upcase(low_memory_mode) = "Y" then close (sigma_e_trace);
		else do;
		
			create (cats(outname, "_sigma_e")) from sigma_e_trace[colname=sigma_e_trace_names];
				append from sigma_e_trace;
			close (cats(outname, "_sigma_e"));
		end;
		
		**MCMC U matrices;
		if upcase(save_mcmc_u) = "Y" then do;
		
			if upcase(low_memory_mode) = "Y" then close (u_matrices_mcmc);
			else do;
			
				u_mcmc_out = repeat(saved_u_traces, j(1, num_mcmc_u, num_subjects))` || u_matrices_mcmc;
				
				create (cats(outname, "_u_mcmc")) from u_mcmc_out[colname=("iteration" || u_mcmc_names)];
					append from u_mcmc_out;
				close (cats(outname, "_u_mcmc"));
				
				free u_mcmc_out;
			end;
		end;
		
		**Post-MCMC conditional U matrices;
		if num_post > 0 then do;
		
			if upcase(low_memory_mode) = "Y" then close (u_matrices_post);
			else do;
			
				u_post_out = repeat(1:num_post, j(1, num_post, num_subjects))` || u_matrices_post;
				
				create (cats(outname, "_u_post")) from u_post_out[colname=("post_mcmc_iteration" || u_post_names)];
					append from u_post_out;
				close (cats(outname, "_u_post"));
				
				free u_post_out;
			end;
		end;
		
		**never-consumer parameters;
		if has_never_consumers = 1 then do;
		
			**alpha1;
			_VARIABLE_ = variable_names[1];
			
			alpha1_mean_t = alpha1_mean`;
			
			create (cats(outname, "_alpha1_mean")) from alpha1_mean_t[rowname=_VARIABLE_ colname=never_consumer_coefficient_names];
				append from alpha1_mean_t[rowname=_VARIABLE_];
			close (cats(outname, "_alpha1_mean"));
			
			alpha1_variance_t = alpha1_variance`;
			
			create (cats(outname, "_alpha1_var")) from alpha1_variance_t[rowname=_VARIABLE colname=never_consumer_coefficient_names];
				append from alpha1_variance_t[rowname=_VARIABLE_];
			close (cats(outname, "_alpha1_var"));
			
			
			
			if upcase(low_memory_mode) = "Y" then close (alpha1_trace);
			else do;
			
				create (cats(outname, "_alpha1")) from alpha1_trace[colname=alpha1_trace_names];
					append from alpha1_trace;
				close (cats(outname, "_alpha1"));
			end;
		
			**consumer probabilities;
			_VARIABLE_ = variable_names[1];
			
			create (cats(outname, "_conprob1_mean")) from consumer_probabilities_mean[rowname=_VARIABLE_ colname={"prob"}];
				append from consumer_probabilities_mean[rowname=_VARIABLE_];
			close (cats(outname, "_conprob1_mean"));
			
			create (cats(outname, "_conprob1_var")) from consumer_probabilities_variance[rowname=_VARIABLE_ colname={"prob"}];
				append from consumer_probabilities_variance[rowname=_VARIABLE_];
			close (cats(outname, "_conprob1_var"));
			
			if upcase(low_memory_mode) = "Y" then close (consumer_probabilities_trace);
			else do;
			
				create (cats(outname, "_conprob1")) from consumer_probabilities_trace[colname=conprob1_trace_name];
					append from consumer_probabilities_trace;
				close (cats(outname, "_conprob1"));
			end;
		end;
		
		if symgetn("DEBUG") then print "end module: mcmc_main_loop";	
	finish;
	
	reset storage=WORK.mcmc_modules;
	
	store module=_all_;
quit;