/* 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.
*/

/*
Program Purpose:
IML modules for updating MCMC parameters using complete conditionals
*/

proc iml;

	**Generates a vector of truncated normal distribution values given a vector of truncation values using rejection sampler (Robert, 1995);
	**Paper source: Robert (1995). Simulation of truncated normal variables. Statistics and Computing, 5, 121-125. https://doi.org/10.48550/arXiv.0907.4010;
	start truncated_normal(trunc_values);
	
	if symgetn("DEBUG") then print "begin module: truncated_normal";		
	
		**number of truncated normal values;
    num_trunc_values = nrow(trunc_values);
      	
    **positive and negative truncation values;
    positive = loc(trunc_values >= 0);
    negative = loc(trunc_values < 0);
      	
    **initialize full truncated normal vector;
    trunc_normal_values = j(num_trunc_values, 1, .);
      	
    **initialize vector of accepted values;
    accepted = j(num_trunc_values, 1, 0);
      	
    **use exponential rejection sampler for any positive truncation values;
    if(^IsEmpty(positive)) then do;
      	
			**subset positive truncation values;
      trunc_values_positive = trunc_values[positive];
      		
      **scale factor for exponential proposal distribution;
      exp_scale_factor = (trunc_values_positive + sqrt(4 + (trunc_values_positive##2)))/2;
      		
      **initialize result vector of truncated normal values for positive truncation values;
      trunc_normal_positive = trunc_normal_values[positive];
      		
      **initialize vector of accepted values for positive truncation values;
      accepted_positive = accepted[positive];
      		
      **exponential rejection sampler for positive truncation values;
      do while(any(^accepted_positive));
      		
       	**determine which values to replace;
      	replace_positive = loc(^accepted_positive);
	      num_to_replace_positive = sum(^accepted_positive);
	      		
	      **calculate candidate values;
	      random_uniform_values = j(num_to_replace_positive, 1, .);
	      call randgen(random_uniform_values, "uniform");
	      candidate_values = trunc_values_positive[replace_positive] - ((1/exp_scale_factor[replace_positive]) # log(random_uniform_values));
	      		
	      **calculate acceptance probabilities for the candidate values and either accept or reject them;
	      acceptance_probability = exp(-0.5 # (candidate_values - exp_scale_factor[replace_positive])##2);
	      call randgen(random_uniform_values, "uniform");
	      accepted_positive[replace_positive] = random_uniform_values < acceptance_probability;
	      		
	      trunc_normal_positive[replace_positive] = candidate_values;
      end;
      		
      trunc_normal_values[positive] = trunc_normal_positive;
    end;
      	
    **use normal rejection sampler for any negative truncation values;
    if(^IsEmpty(negative)) then do;
      	
      **subset negative truncation values;
      trunc_values_negative = trunc_values[negative];
      		
      **initialize result vector of truncated normal values for negative truncation values;
      trunc_normal_negative = trunc_normal_values[negative];
      		
      **initialize vector of accepted values for positive truncation values;
      accepted_negative = accepted[negative];
      		
      **normal rejection sampler for negative truncation values;
      do while(any(^accepted_negative));
      		
      	**determine which values to replace;
      	replace_negative = loc(^accepted_negative);
	      num_to_replace_negative = sum(^accepted_negative);
	      		
	      **draw from a normal distribution with a mean of 0 and a standard deviation of 1;
	      new_trunc_normal_negative = j(num_to_replace_negative, 1, .);
	      call randgen(new_trunc_normal_negative, "normal");
	      accepted_negative[replace_negative] = new_trunc_normal_negative >= trunc_values_negative[replace_negative];
	      		
	      trunc_normal_negative[replace_negative] = new_trunc_normal_negative;    		
      end;
      		
      trunc_normal_values[negative] = trunc_normal_negative;
    end;
      	
	if symgetn("DEBUG") then print "end module: truncated_normal";
    return(trunc_normal_values);
	finish;
	
	**Updates W matrix for episodically consumed foods;
	**Appendix A.9 and A.10 in Zhang, et al. (2011);
	start update_w_matrix(w_matrix,
												xbeta_u,
												sigma_e,
												recall_availability,
												episodic_indicator_matrices,
												num_subjects,
												num_episodic,
												num_variables,
												num_recalls);
		
		if symgetn("DEBUG") then print "begin module: update_w_matrix";
		inv_sigma_e = inv(sigma_e);
		
		w_matrix_updated = w_matrix;
		do day_k = 1 to num_recalls;
		
			w_matrix_updated_k = w_matrix_updated$day_k;
			episodic_indicator_matrix_k = episodic_indicator_matrices$day_k;
			xbeta_u_k = xbeta_u$day_k;
			
			do var_j = 1 to 2*num_episodic;
			
				food_j = int((var_j + 1)/2);
				observed = episodic_indicator_matrix_k[,food_j];
				
				**C1 and C2 in W matrix complete conditionals from Appendix A.9 and A.10 of Zhang, et al. (2011);
				c1 = j(num_subjects, 1, 0);
				do var_jj = 1 to num_variables;
				
					if var_jj = var_j then do;
					
						c1 = c1 + xbeta_u_k[,var_jj] # inv_sigma_e[var_jj, var_jj];
					end;
					else do;
					
						c1 = c1 - (w_matrix_updated_k[,var_jj] - xbeta_u_k[,var_jj]) # inv_sigma_e[var_j, var_jj];
					end;
				end;
				
				c2 = 1/inv_sigma_e[var_j, var_j];
				
				mu = c2 # c1;
				sigma = sqrt(c2);
				
				if mod(var_j, 2) = 1 then do;
				
					**Update W matrix column for indicator using the complete conditional from Appendix A.9 of Zhang, et al. (2011);
					truncation_values = (2 # observed - 1) # -mu/sigma;
					trunc_normal = truncated_normal(truncation_values);
					
					if any(observed) then do;
					
						w_matrix_updated_k[loc(observed), var_j] = (mu[loc(observed)] + sigma # trunc_normal[loc(observed)]) # recall_availability[loc(observed), day_k];
					end;
					
					if any(^observed) then do;
					
						w_matrix_updated_k[loc(^observed), var_j] = (mu[loc(^observed)] - sigma # trunc_normal[loc(^observed)]) # recall_availability[loc(^observed), day_k];
					end;
				end;
				else do;
				
					**Update the W matrix column for amount using the complete conditional from Appendix A.10 of Zhang, et al. (2011);
					normals = j(num_subjects, 1, .);
					call randgen(normals, "Normal");
					normals = mu + normals # sigma;
					
					if any(^observed) then do;
					
						w_matrix_updated_k[loc(^observed), var_j] = normals[loc(^observed)] # recall_availability[loc(^observed), day_k];
					end;
				end;
			end;
			
			w_matrix_updated$day_k = w_matrix_updated_k;
		end;
		
		if symgetn("DEBUG") then print "end module: update_w_matrix";
		return w_matrix_updated;
	finish;
	
	**Weighted W-XBeta-U multiplied by unweighted W-XBeta-U for all variable combinations, summed over all recalls;
	**W-XBeta-U corresponds to the error (epsilon) term in Equation 3.5 of Zhang, et al. (2011);
	**Pre-calculated component in complete conditionals for r, theta, and V elements;
	**See Appendix A.5 of Zhang, et al. (2011);
	**This is calculated before updating r, theta, and V in order to save time - sigma-e varies with each element, but the error terms change just once per iteration;
	start calculate_w_cross_residual_sum(w_matrix,
																			 xbeta_u,
																			 recall_availability,
																			 subject_weighting,
																			 num_variables,
																			 num_recalls);
																			 
		if symgetn("DEBUG") then print "begin module: calculate_w_cross_residual_sum";		
        alg = symgetn("v143code") + 1;
		
		w_cross_residual_sum = j(num_variables, num_variables, 0);
		
		do day_k = 1 to num_recalls;
		
			**calculate residual error terms W - XBeta - U;
			w_minus_xbeta_u = (w_matrix$day_k - xbeta_u$day_k) # recall_availability[,day_k];
			
			**sum the cross products of the residuals;
/* BW CHANGE 2025-10-16 */
            tmp = (w_minus_xbeta_u # subject_weighting)`;
			w_cross_residual_sum = w_cross_residual_sum + matMultAB(tmp, w_minus_xbeta_u, alg);
			*w_cross_residual_sum = w_cross_residual_sum + (w_minus_xbeta_u # subject_weighting)` * w_minus_xbeta_u;
					
		end;
		
		if symgetn("DEBUG") then print "end module: calculate_w_cross_residual_sum";
		return w_cross_residual_sum;
	finish;
	
	**Updates r matrix using Metropolis-Hastings step;
	**See Appendix A.5 of Zhang, et al. (2011);
	start update_r_matrix(r_matrix,
												theta_matrix,
												v_matrix,
												w_cross_residual_sum,
												recall_availability,
												subject_weighting);
												
		if symgetn("DEBUG") then print "begin module: update_r_matrix";
		alg = symgetn("v143code") + 1;
		
		**if r matrix is empty, then return it immediately (this is because there aren't two or more episodic foods);
		if IsEmpty(r_matrix) then do;
		
			return r_matrix;
		end;
		
		**uniform spacing between possible r values;
		r_spacing = 2*0.99/200;
		
		**discrete list of possible r-values (see Appendix A.5 of Zhang, et al. (2011));
		r_possible = do(-0.99, 0.99, r_spacing);
		
		**initialize current r matrix;
		r_current = r_matrix;
		
		**loop through current r values to update them;
		do i = 1 to nrow(r_current);
		
			**calculate minimum and maximum r;
			r_min = min(r_possible, r_current[i]);
			r_max = max(r_possible, r_current[i]);
			
			**choose candidate r value;
			**1. If r[i] <= r.min, choose candidate from current r[i] and r[i] + 1 or 2 r spacings;
			**2. If r[i] is between r.min and r.max, choose candidate from current r[i] and r[i] +/- 1 r spacing;
			**3. If r[i] >= r.max, choose candidate from current r[i] and r[i] - 1 or 2 r spacings;
			
			**random uniform number to randomly select candidate from options;
			candidate_choice = j(1, 1, .);
			call randgen(candidate_choice, "Uniform");
		
			if r_current[i] <= r_min then do;
			
				candidate = (r_current[i])*(candidate_choice <= 1/3) +
										(r_current[i] + r_spacing)*(candidate_choice > 1/3 & candidate_choice <= 2/3) +
										(r_current[i] + 2*r_spacing)*(candidate_choice > 2/3);
			end;
			else if r_current[i] >= r_max then do;
			
				candidate = (r_current[i])*(candidate_choice <= 1/3) +
										(r_current[i] - r_spacing)*(candidate_choice > 1/3 & candidate_choice <= 2/3) +
										(r_current[i] - 2*r_spacing)*(candidate_choice > 2/3);
			end;
			else do;
			
				candidate = (r_current[i])*(candidate_choice <= 1/3) +
										(r_current[i] + r_spacing)*(candidate_choice > 1/3 & candidate_choice <= 2/3) +
										(r_current[i] - r_spacing)*(candidate_choice > 2/3);
			end;
			
			**if r candidate is between -1 and 1, perform Metropolis-Hastings step to determine whether to accept it;
			**otherwise, automatically reject the candidate and the current r is kept the same;
			if candidate >= -1 & candidate <= 1 then do;
			
				**create candidate r matrix;
				r_candidate = r_current;
				r_candidate[i] = candidate;
				
				**V matrix based on current and candidate r values;

				v_current = update_v_indicator(v_matrix, 
																			 r_current,
																			 theta_matrix);

				v_candidate = update_v_indicator(v_matrix,
																				 r_candidate,
																				 theta_matrix);
																		 
				**sigma-e based on current and candidate r values;
/* BW CHANGE 2025-10-16 */
				*sigma_e_current = v_current * v_current`;	
				*sigma_e_candidate = v_candidate * v_candidate`;
				sigma_e_current   = matMultAAt(v_current, alg);	
				sigma_e_candidate = matMultAAt(v_candidate, alg);
				
				**inverse sigma.e for current and candidate values;
				inv_sigma_e_current = inv(sigma_e_current);
				inv_sigma_e_candidate = inv(sigma_e_candidate);
				
				**weighted sum of number of recalls across all subjects;
				num_recalls_per_subject = recall_availability[,+];
				weighted_sum_recalls = sum(num_recalls_per_subject # subject_weighting);
				
				**current and candidate log-likelihood for Metropolis-Hastings step (Appendix A.5 of Zhang, et al. (2011));
				current_log_likelihood = -0.5*weighted_sum_recalls*log(1 - r_current[i]**2) - 0.5*sum(inv_sigma_e_current # w_cross_residual_sum);
				candidate_log_likelihood = -0.5*weighted_sum_recalls*log(1 - r_candidate[i]**2) - 0.5*sum(inv_sigma_e_candidate # w_cross_residual_sum);
				
				**calculate Metropolis-Hastings acceptance probability;
				if candidate_log_likelihood >= current_log_likelihood then do;
    		
    			acceptance_probability = 1;
    		end;
    		else do;
    		
    			acceptance_probability = exp(candidate_log_likelihood - current_log_likelihood);
    		end;
			
				**current r value is replaced with the candidate with probability equal to the acceptance probability, otherwise the current r value is kept;
				acceptance_choice = j(1, 1, .);
				call randgen(acceptance_choice, "Uniform");
				r_current[i] = r_candidate[i]*(acceptance_choice <= acceptance_probability) + r_current[i]*(acceptance_choice > acceptance_probability);
			end;
		end;
	
		if symgetn("DEBUG") then print "end module: update_r_matrix";
		return r_current;
	finish;
	
	**Updates theta matrix using Metropolis-Hastings step;
	**See Appendix A.5 of Zhang, et al. (2011);
	start update_theta_matrix(theta_matrix,
														r_matrix,
														v_matrix,
														w_cross_residual_sum);
													
		if symgetn("DEBUG") then print "begin module: update_theta_matrix";
		alg = symgetn("v143code") + 1;
		
		**if theta matrix is empty, then return it immediately (this is because there aren't two or more episodic foods);
		if IsEmpty(theta_matrix) then do;
		
			return theta_matrix;
		end;
		
		**uniform spacing between possible theta values;
		pi = arcos(-1);
		theta_spacing = 2*pi*0.99/200;
		
		**discrete list of possible theta values (see Appendix A.5 of Zhang, et al. (2011));
		theta_possible = do(-pi, pi, theta_spacing);
		
		**initialize current theta matrix;
		theta_current = theta_matrix;
		
		**loop through current theta values to update them;
		do i = 1 to nrow(theta_current);
		
			**calculate minimum and maximum theta;
			theta_min = min(theta_possible, theta_current[i]);
			theta_max = max(theta_possible, theta_current[i]);
			
			**choose candidate theta value;
			**1. If theta[i] <= theta.min, choose candidate from current theta[i] and theta[i] + 1 or 2 theta spacings;
			**2. If theta[i] is between theta.min and theta.max, choose candidate from current theta[i] and theta[i] +/- 1 theta spacing;
			**3. If theta[i] >= theta.max, choose candidate from current theta[i] and theta[i] - 1 or 2 theta spacings;
			
			**random uniform number to randomly select candidate from options;
			candidate_choice = j(1, 1, .);
			call randgen(candidate_choice, "Uniform");
			
			if theta_current[i] <= theta_min then do;
										
				candidate = (theta_current[i])*(candidate_choice <= 1/3) +
										(theta_current[i] + theta_spacing)*(candidate_choice > 1/3 & candidate_choice <= 2/3) +
										(theta_current[i] + 2*theta_spacing)*(candidate_choice > 2/3);
			end;
			else if theta_current[i] >= theta_max then do;
			
				candidate = (theta_current[i])*(candidate_choice <= 1/3) +
										(theta_current[i] - theta_spacing)*(candidate_choice > 1/3 & candidate_choice <= 2/3) +
										(theta_current[i] - 2*theta_spacing)*(candidate_choice > 2/3);
			end;
			else do;
			
				candidate = (theta_current[i])*(candidate_choice <= 1/3) +
										(theta_current[i] + theta_spacing)*(candidate_choice > 1/3 & candidate_choice <= 2/3) +
										(theta_current[i] - theta_spacing)*(candidate_choice > 2/3);
			end;
			
			**create candidate theta matrix;
			theta_candidate = theta_current;
			theta_candidate[i] = candidate;
			
			**V matrix based on current and candidate theta values;
			v_current = update_v_indicator(v_matrix,
																		 r_matrix,
																		 theta_current);
			v_candidate = update_v_indicator(v_matrix,
																			 r_matrix,
																			 theta_candidate);
																			 
			**sigma-e based on current and candidate theta values;
/* BW CHANGE 2025-10-16 */
			*sigma_e_current = v_current * v_current`;
			*sigma_e_candidate = v_candidate * v_candidate`;
			sigma_e_current   = matMultAAt(v_current, alg);
			sigma_e_candidate = matMultAAt(v_candidate, alg);
			
			**inverse sigma-e for current and candidate values;
			inv_sigma_e_current = inv(sigma_e_current);
			inv_sigma_e_candidate = inv(sigma_e_candidate);
			
			**current and candidate log-likelihood for Metropolis-Hastings step (Appendix A.5 of Zhang, et al. (2011));
			current_log_likelihood = -0.5*sum(inv_sigma_e_current # w_cross_residual_sum);
			candidate_log_likelihood = -0.5*sum(inv_sigma_e_candidate # w_cross_residual_sum);
			
			**calculate Metropolis-Hastings acceptance probability;
			if candidate_log_likelihood >= current_log_likelihood then do;
    		
    		acceptance_probability = 1;
    	end;
    	else do;
    		
    		acceptance_probability = exp(candidate_log_likelihood - current_log_likelihood);
    	end;
			
			**current theta value is replaced with the candidate with probability equal to the acceptance probability, otherwise the current theta value is kept;
			acceptance_choice = j(1, 1, .);
			call randgen(acceptance_choice, "Uniform");
			theta_current[i] = theta_candidate[i]*(acceptance_choice <= acceptance_probability) + theta_current[i]*(acceptance_choice > acceptance_probability);
		end;
		
		if symgetn("DEBUG") then print "end module: update_theta_matrix";
		return theta_current;
	finish;
	
	**Updates V matrix in three parts;
	**1. Diagonal elements for amount variables;
	**2. Off-diagonal elements for amount variables;
	**3. Elements for indicator variables;
	start update_v_matrix(v_matrix,
												r_matrix,
												theta_matrix,
												w_cross_residual_sum,
												recall_availability,
												subject_weighting,
												num_episodic,
												num_daily);
		
		if symgetn("DEBUG") then print "begin module: update_v_matrix";
		v_updated = v_matrix;
	
		**1. Diagonal elements for amount variables;
		v_updated = update_v_amount_diagonal(v_updated,
																				 r_matrix,
																				 theta_matrix,
																				 w_cross_residual_sum,
																				 recall_availability,
																				 subject_weighting,
																				 num_episodic,
																				 num_daily);
		
		**2. Off-diagonal elements for amount variables;
		v_updated = update_v_amount_off_diagonal(v_updated,
																						 r_matrix,
																						 theta_matrix,
																						 w_cross_residual_sum,
																						 num_episodic,
																						 num_daily);
		
		**3. Elements for indicator variables;
		v_updated = update_v_indicator(v_updated,
																	 r_matrix,
																	 theta_matrix);
		
		if symgetn("DEBUG") then print "end module: update_v_matrix";
		return v_updated;
	finish;
	
	**Updates V matrix diagonal elements for amount variables using Metropolis-Hastings step;
	**See Vqq complete conditional from Zhang, et al. (2011);
	start update_v_amount_diagonal(v_matrix,
																 r_matrix,
																 theta_matrix,
																 w_cross_residual_sum,
																 recall_availability,
																 subject_weighting,
																 num_episodic,
																 num_daily);
		
		if symgetn("DEBUG") then print "begin module: update_v_amount_diagonal";
		alg = symgetn("v143code") + 1;
		
		**initialize current V matrix;
		v_current = v_matrix;
		
		**loop through amount diagonals and update them;
		do var_num = 1 to num_episodic + num_daily;
		
			**index of amount variable;
			if var_num <= num_episodic then do;
			
				i = 2*var_num;
			end;
			else do;
			
				i = num_episodic + var_num;
			end;
			
			**extract candidate value for V[i,i];
			current = v_current[i,i];
			candidate = j(1, 1, .);
			call randgen(candidate, "Uniform", current-0.2, current+0.2);
			
			**if candidate is between -3 and 3, perform Metropolis-Hastings step to determine whether to accept it;
    	**otherwise, automatically reject the candidate and the current value is kept;
    	if candidate >= -3 & candidate <= 3 then do;
    	
    		**replace V[i,i] in candidate V matrix with candidate value;
    		v_candidate = v_current;
    		v_candidate[i,i] = candidate;
    		
    		**proposed v matrices based on current and candidate values;
    		v_proposed_current = update_v_indicator(v_current,
    																						r_matrix,
    																						theta_matrix);
    		v_proposed_candidate = update_v_indicator(v_candidate,
    																							r_matrix,
    																							theta_matrix);
    																						
    		**proposed sigma-e from current and candidate values;
/* BW CHANGE 2025-10-16 */
    		*sigma_e_current = v_proposed_current * v_proposed_current`;
    		*sigma_e_candidate = v_proposed_candidate * v_proposed_candidate`;
			sigma_e_current   = matMultAAt(v_proposed_current, alg);
    		sigma_e_candidate = matMultAAt(v_proposed_candidate, alg);
    		
    		**inverse of the proposed sigma-e;
    		inv_sigma_e_current = inv(sigma_e_current);
    		inv_sigma_e_candidate = inv(sigma_e_candidate);
    		
    		**weighted sum of number of recalls across all subjects;
    		num_recalls_per_subject = recall_availability[,+];
    		weighted_sum_recalls = sum(num_recalls_per_subject # subject_weighting);
    		
    		**current and candidate log-likelihood for Metropolis-Hastings step (Appendix A.5 of Zhang, et al. (2011));
    		current_log_likelihood = -0.5*weighted_sum_recalls*log(v_current[i,i]**2) - 0.5*sum(inv_sigma_e_current # w_cross_residual_sum);
    		candidate_log_likelihood = -0.5*weighted_sum_recalls*log(v_candidate[i,i]**2) - 0.5*sum(inv_sigma_e_candidate # w_cross_residual_sum);
    		
    		**calculate Metropolis-Hastings acceptance probability;
    		if candidate_log_likelihood >= current_log_likelihood then do;
    		
    			acceptance_probability = 1;
    		end;
    		else do;
    		
    			acceptance_probability = exp(candidate_log_likelihood - current_log_likelihood);
    		end;
				
				**current value is replaced with the candidate with probability equal to the acceptance probability, otherwise the current value is kept;
				acceptance_choice = j(1, 1, .);
				call randgen(acceptance_choice, "Uniform");
				v_current[i,i] = v_candidate[i,i]*(acceptance_choice <= acceptance_probability) + v_current[i,i]*(acceptance_choice > acceptance_probability);
    	end;
		end;
		
		if symgetn("DEBUG") then print "end module: update_v_amount_diagonal";
		return v_current;
	finish;
	
	**Updates V matrix off-diagonal elements for amount variables using Metropolis-Hastings step;
	**see sample Vpq complete conditional from Zhang, et al. (2011);
	start update_v_amount_off_diagonal(v_matrix,
																		 r_matrix,
																		 theta_matrix,
																		 w_cross_residual_sum,
																		 num_episodic,
																		 num_daily);
			
		if symgetn("DEBUG") then print "begin module: update_v_amount_off_diagonal";
		alg = symgetn("v143code") + 1;
		
		**initialize current V matrix;
		v_current = v_matrix;
		
		**loop through every row of V matrix corresponding to amount variables and update off-diagonal elements;
		do var_num = 1 to num_episodic + num_daily;
		
			**index of amount variable and number of elements to update;
			if var_num <= num_episodic then do;
			
				i = 2*var_num;
				num_off_diag = i - 2;
			end;
			else do;
			
				i = num_episodic + var_num;
				num_off_diag = i - 1;
			end;
		
			**loop through off-diagonal elements and update them;
			do j = 1 to num_off_diag;
			
				**Extract candidate value for V[i,j];
				current = v_current[i,j];
				candidate = j(1, 1, .);
				call randgen(candidate, "Uniform", current-0.2, current+0.2);
				
				**if candidate is between -3 and 3, perform Metropolis-Hastings step to determine whether to accept it;
      	**otherwise, automatically reject the candidate and the current value is kept;
      	if candidate >= -3 & candidate <= 3 then do;
      	
      		**replace V[i,j] in candidate V matrix with candidate value;
      		v_candidate = v_current;
      		v_candidate[i,j] = candidate;
      		
      		**proposed v matrices based on current and candidate values;
      		v_proposed_current = update_v_indicator(v_current,
    																							r_matrix,
    																							theta_matrix);
    			v_proposed_candidate = update_v_indicator(v_candidate,
    																								r_matrix,
    																								theta_matrix);
    																							
    			**proposed sigma-e from current and candidate values;
/* BW CHANGE 2025-10-16 */
    			*sigma_e_current = v_proposed_current * v_proposed_current`;
    			*sigma_e_candidate = v_proposed_candidate * v_proposed_candidate`;
				sigma_e_current   = matMultAAt(v_proposed_current, alg);
    			sigma_e_candidate = matMultAAt(v_proposed_candidate, alg);
    		
    			**inverse of the proposed sigma-e;
    			inv_sigma_e_current = inv(sigma_e_current);
    			inv_sigma_e_candidate = inv(sigma_e_candidate);
    			
    			**current and candidate log-likelihood for Metropolis-Hastings step (Appendix A.5 of Zhang, et al. (2011));
    			current_log_likelihood = -0.5*sum(inv_sigma_e_current # w_cross_residual_sum);
    			candidate_log_likelihood = -0.5*sum(inv_sigma_e_candidate # w_cross_residual_sum);
    			
    			**calculate Metropolis-Hastings acceptance probability;
					if candidate_log_likelihood >= current_log_likelihood then do;
    		
    				acceptance_probability = 1;
    			end;
    			else do;
    		
    				acceptance_probability = exp(candidate_log_likelihood - current_log_likelihood);
    			end;
				
					**current value is replaced with the candidate with probability equal to the acceptance probability, otherwise the current value is kept;
					acceptance_choice = j(1, 1, .);
					call randgen(acceptance_choice, "Uniform");
					v_current[i,j] = v_candidate[i,j]*(acceptance_choice <= acceptance_probability) + v_current[i,j]*(acceptance_choice > acceptance_probability);
      	end;
			end;
		end;
		
		if symgetn("DEBUG") then print "end module: update_v_amount_off_diagonal";
		return v_current;
	finish;
	
	**Updating V matrix indicator variables for both the main Gibbs sampler and for the current/candidate sigma-e in the Metropolis-Hastings steps for r, theta, and v elements;
	**See section 3.2 of Zhang, et al. (2011) for the formula used in this function;
	start update_v_indicator(v_matrix,
													 r_matrix,
													 theta_matrix);
						
		if symgetn("DEBUG") then print "begin module: update_v_indicator";
		**updates are only needed if there are 2 or more episodic variables;
		if IsEmpty(r_matrix) then do;
		
			return v_matrix;
		end;
		
		**initialize update V matrix to current V matrix;
		v_updated = v_matrix;
		
		do q = 1 to nrow(r_matrix);
		
			do p = 1 to 2*q - 1;
			
				if p = 1 then do;
				
					v_updated[2*q + 1, p] = r_matrix[q] * sin(theta_matrix[p + (q-1)**2]);
				end;
				else do;
				
					v_updated[2*q + 1, p] = r_matrix[q] * prod(cos(theta_matrix[1:(p-1) + (q-1)**2])) * sin(theta_matrix[p + (q-1)**2]);
				end;
				
				v_updated[2*q + 1, 2*q] = r_matrix[q] * prod(cos(theta_matrix[1:p + (q-1)**2]));
			end;
			
			v_updated[2*q + 1, 2*q + 1] = sqrt(1 - r_matrix[q]**2);
			v_updated[2*q + 2, 2*q + 1] = -sum(v_updated[2*q + 1, 1:(2*q)] # v_updated[2*q + 2, 1:(2*q)])/v_updated[2*q + 1, 2*q + 1];
		end;
		
		if symgetn("DEBUG") then print "end module: update_v_indicator";
		return v_updated;
	finish;
	
	**Updates Sigma-u matrix;
	**Uses complete conditional in Appendix A.6 of Zhang, et al. (2011);
	start update_sigma_u(sigma_u,
											 sigma_u_prior,
											 u_matrix,
											 subject_weighting,
											 num_subjects,
											 num_variables,
											 sigma_u_constant);
					
		if symgetn("DEBUG") then print "begin module: update_sigma_u";
		alg = symgetn("v143code") + 1;
		
		**check here for constant sigma-u parameter - return sigma-u unchanged if it is true;
		if upcase(sigma_u_constant) = "Y" then do;
		
			return sigma_u;
		end;
	
		**sigma-u prior degrees of freedom;
		sigma_u_prior_df = nrow(sigma_u_prior) + 2;
		
		**calculate inverse Wishart matrix parameter;
/* BW CHANGE 2025-10-16 */
		*inv_wishart = (sigma_u_prior_df - nrow(sigma_u) - 1) # sigma_u_prior + (subject_weighting # u_matrix)` * u_matrix;
		tmp         = (subject_weighting # u_matrix)`;
        tmp         = matMultAB(tmp, u_matrix, alg);
		inv_wishart = (sigma_u_prior_df - nrow(sigma_u) - 1) # sigma_u_prior + tmp;
	
		**calculate inverse Wishart degrees of freedom;
		wishart_df = sigma_u_prior_df + num_subjects;
		
		**Adjust inverse Wishart matrix diagonals and degrees of freedom by 0.0011 times the number of subjects;
		inv_wishart = inv_wishart + 0.0011*num_subjects # i(ncol(sigma_u));
		wishart_df = wishart_df + 0.0011*num_subjects;
		
		**convert inverse Wishart matrix into a symmetric matrix;
		inv_wishart = (inv_wishart + inv_wishart`)/2;
	
		**take the inverse of the inverse Wishart matrix to compute the matrix parameter for the forward Wishart function;
		fwd_wishart = inv(inv_wishart);

		**calculate inverse sigma-u using the forward Wishart function;
/* BW CHANGE 2025-10-17 */
		*inv_sigma_u = randwishart(1, wishart_df, fwd_wishart);
		*inv_sigma_u = shape(inv_sigma_u, num_variables, num_variables);
		inv_sigma_u = myRandWishart(wishart_df, fwd_wishart, alg);

		**calculate sigma-u;
		sigma_u_updated = inv(inv_sigma_u);
		
		if symgetn("DEBUG") then print "end module: update_sigma_u";
		return sigma_u_updated;
	finish;
	
	**Updates U matrix;
	**Uses complete conditional in Appendix A.8 of Zhang, et al. (2011);
	**Uses a Metropolis-Hastings step to sample values for never-consumers when they are present;
	start 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 symgetn("DEBUG") then print "begin module: update_u_matrix";
		alg = symgetn("v143code") + 1;
		
		**initialize candidate U matrix to input U matrix;
		u_matrix_current = u_matrix;
		u_matrix_candidate = u_matrix;
		
		**calculating inverse of sigma-e and sigma-u;
		inv_sigma_e = inv(sigma_e);
		inv_sigma_u = inv(sigma_u);
  	
		**calculate total recalls per subject;
		recalls_per_subject = recall_availability[,+];
	
		**Calculating C1 in complete conditional;
		w_minus_xbeta_sum = j(num_subjects, num_variables, 0);
		do day_k = 1 to num_recalls;
		
			w_minus_xbeta_sum = w_minus_xbeta_sum + (w_matrix$day_k - xbeta$day_k) # recall_availability[,day_k];
		end;

/* BW CHANGE: 2025-10-17 */
		c1 = matMultAB(w_minus_xbeta_sum, inv_sigma_e, alg);
		*c1 = w_minus_xbeta_sum * inv_sigma_e;
	
		**calculate candidate U matrix for each number of recalls;
		do count_k = 1 to num_recalls;
		
			if any(recalls_per_subject = count_k) then do;
			
				**calculating C2 in complete conditional;
				c2 = inv(inv_sigma_u + count_k # inv_sigma_e);
				
				**calculate the mean of the candidate U matrix;
/* BW CHANGE 2025-10-17 */
				u_matrix_mean = matMultAB(c1[loc(recalls_per_subject = count_k),], c2, alg);
				*u_matrix_mean = c1[loc(recalls_per_subject = count_k),] * c2;
				
				**calculate standard deviation of candidate U matrix by taking the matrix square root of C2;
/* BW CHANGE 2025-10-17 */
				*call eigen(c2_eigvals, c2_eigvecs, c2);
				tmp = eigenDecomp(c2, alg);
				c2_eigvals = tmp$"eigenvalues";
				c2_eigvecs = tmp$"eigenvectors";

/* BW CHANGE 2025-10-17 */
				*u_matrix_std_dev = c2_eigvecs * diag(sqrt(c2_eigvals)) * c2_eigvecs`;
				tmp = sqrt(c2_eigvals);
				u_matrix_std_dev = matMultADAt(c2_eigvecs, tmp, alg);
				free tmp;
				
				**update the candidate U matrix for the current number of recalls according to the complete conditional in Appendix A.8 of Zhang, et al. (2011);
				**normal distribution with a mean of C1 * C2 and a standard deviation of sqrt(C2);
				normals = j(sum(recalls_per_subject = count_k), num_variables, .);
				call randgen(normals, "Normal");
		
/* BW CHANGE 2025-10-17 */
				*u_matrix_candidate[loc(recalls_per_subject = count_k),] = u_matrix_mean + normals * u_matrix_std_dev;
				u_matrix_candidate[loc(recalls_per_subject = count_k),] = u_matrix_mean + matMultAB(normals, u_matrix_std_dev, alg);

			end;
		end;
	
		**if never-consumers are allowed, perform Metropolis-Hastings step for the first episodic variable for each never-consumer and automatically accept consumers;
		**otherwise, accept entire candidate U matrix as-is;
		if has_never_consumers = 1 then do;
		
			if any(conni1 <= 0) then do;
			
				**never-consumer subjects;
				never_consumers = (conni1 <= 0);
				num_never_consumers = sum(never_consumers);
				
				**calculate log-likelihood for current and candidate U matrices;
				log_likelihood_candidate = j(num_subjects, 1, 0);
				log_likelihood_current = j(num_subjects, 1, 0);
				
				do day_k = 1 to num_recalls;
				
					xbeta_k = xbeta$day_k;
					
					xbeta_u_candidate = xbeta_k[loc(never_consumers),1] + u_matrix_candidate[loc(never_consumers),1];
					xbeta_u_current = xbeta_k[loc(never_consumers),1] + u_matrix_current[loc(never_consumers),1];
					
					density_candidate = (1 - cdf("Normal", xbeta_u_candidate)) <> 1e-8;
					density_current = (1 - cdf("Normal", xbeta_u_current)) <> 1e-8;
					
					log_likelihood_candidate[loc(never_consumers)] = log_likelihood_candidate[loc(never_consumers)] - recall_availability[loc(never_consumers),day_k] # log(density_candidate);
					log_likelihood_current[loc(never_consumers)] = log_likelihood_current[loc(never_consumers)] - recall_availability[loc(never_consumers),day_k] # log(density_current);
				end;
				
				**calculate acceptance probabilites for each subject;
				**log-likelihood ratios above zero are set to zero (corresponding to an acceptance probability of 1);
				log_likelihood_ratio = (log_likelihood_candidate - log_likelihood_current) >< 0;
				acceptance_probability = exp(log_likelihood_ratio);
				
				**accept the candidate for each subject with a probability equal to the acceptance probability;
				**otherwise reject the candidate and use the original U matrix value;
				acceptance_choice = j(num_subjects, 1, .);
				call randgen(acceptance_choice, "Uniform");
				u_matrix_updated = u_matrix_candidate # (acceptance_choice <= acceptance_probability) + u_matrix_current # (acceptance_choice > acceptance_probability);
			end;
			else do;
			
				u_matrix_updated = u_matrix_candidate;
			end;
		end;
		else do;
		
			u_matrix_updated = u_matrix_candidate;
		end;
	
		if symgetn("DEBUG") then print "end module: update_u_matrix";
		return u_matrix_updated;
	finish;
	
	**Updates betas for each variable;
	**Uses complete conditional defined in Appendix A.7 of Zhang, et al. (2011);
	**If never-consumers are allowed, uses a Metropolis-Hastings step to sample beta values for the first episodic variable;
	start 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,
										beta1);
						
		if symgetn("DEBUG") then print "begin module: update_beta";
		alg = symgetn("v143code") + 1;
		
		inv_sigma_e = inv(sigma_e);
		
		**update betas for each episodic and daily variable;
		beta_updated = ListCreate(num_variables);
		
		do var_j = 1 to num_variables;
	
			**number of covariates for the current variable;
			num_coefficients_j = num_coefficients[var_j];
			
			**calculate the C1 term in the beta complete conditional from Appendix A.7 of Zhang, et al. (2011);
/* BW CHANGE 2025-10-22 */
			*c1 = inv(beta_covariance_prior$var_j) * beta_mean_prior$var_j;
			c1 = matMultAB(inv(beta_covariance_prior$var_j), beta_mean_prior$var_j, alg);
		
			do var_jj = 1 to num_variables;
			
				do day_k = 1 to num_recalls;
				
					w_matrix_k = w_matrix$day_k;
					xbeta_k = xbeta$day_k;
				
					if var_jj = var_j then do;
					
						c1 = c1 + (wt_covariate_matrices$var_j$day_k * (w_matrix_k[,var_j] - u_matrix[,var_j])) # inv_sigma_e[var_j,var_j];
					end;
					else do;
					
						c1 = c1 + (wt_covariate_matrices$var_j$day_k * (w_matrix_k[,var_jj] - xbeta_k[,var_jj] - u_matrix[,var_jj])) # inv_sigma_e[var_jj,var_j];
					end;
				end;
			end;
			
			**calculate the C2 term in the beta complete conditional from Appendix A.7 of Zhang, et al. (2011);
			inverse_c2 = inv(beta_covariance_prior$var_j) + inv_sigma_e[var_j,var_j] # wt_covariate_sq_sums$var_j;
			c2 = inv(inverse_c2);
			
			**if never-consumers are allowed, perform Metropolis-Hastings step to update beta for first episodic variable;
			**otherwise, use the complete conditional in Appendix A.7 of Zhang, et al. (2011);
			if has_never_consumers = 1 & var_j = 1 then do;
			
				**never-consumer subjects;
				never_consumers = (conni1 <= 0);
				
				**current beta for first episodic variable;
				beta1_current = beta1;
				
				**candidate beta for first episodic variable;
				beta1_variance = 2 # c2;
/* BW CHANGE 2025-10-22 */
				*call eigen(beta1_variance_eigvals, beta1_variance_eigvecs, beta1_variance);
				tmp = eigenDecomp(beta1_variance, alg);
				beta1_variance_eigvals = tmp$"eigenvalues";
				beta1_variance_eigvecs = tmp$"eigenvectors";

				*beta1_std_dev = beta1_variance_eigvecs * diag(sqrt(beta1_variance_eigvals)) * beta1_variance_eigvecs`;
				tmp = sqrt(beta1_variance_eigvals);
				beta1_std_dev = matMultADAt(beta1_variance_eigvecs, tmp, alg);
				free tmp;
				
				normals = j(num_coefficients_j, 1, .);
				call randgen(normals, "Normal");
				
/* BW CHANGE 2025-10-22 */
				*beta1_candidate = beta1_current + beta1_std_dev * normals;
				beta1_candidate = beta1_current + matMultAB(beta1_std_dev, normals, alg);
				
				**calculate log-likelihood for candidate and current beta;
/* BW CHANGE 2025-10-22 */
				*log_likelihood_candidate = c1` * beta1_candidate - 0.5*(beta1_candidate` * inverse_c2 * beta1_candidate);
				*log_likelihood_current = c1` * beta1_current - 0.5*(beta1_current` * inverse_c2 * beta1_current);
				log_likelihood_candidate = matMultAB(c1`, beta1_candidate, alg) - 0.5*matMultABAt(beta1_candidate`, inverse_c2, alg);
				log_likelihood_current   = matMultAB(c1`, beta1_current, alg)   - 0.5*matMultABAt(beta1_current`,   inverse_c2, alg);
				
				if any(never_consumers) then do;
				
					do day_k = 1 to num_recalls;
					
						covariate_matrix1_k = covariate_matrices$1$day_k;
						
						xbeta_u1_candidate = covariate_matrix1_k[loc(never_consumers),] * beta1_candidate + u_matrix[loc(never_consumers),1];
						xbeta_u1_current = covariate_matrix1_k[loc(never_consumers),] * beta1_current + u_matrix[loc(never_consumers),1];
						
						density_candidate = (1 - cdf("Normal", xbeta_u1_candidate)) <> 1e-8;
						density_current = (1 - cdf("Normal", xbeta_u1_current)) <> 1e-8;
						
						log_likelihood_candidate = log_likelihood_candidate - sum(recall_availability[loc(never_consumers),day_k] # log(density_candidate));
						log_likelihood_current = log_likelihood_current - sum(recall_availability[loc(never_consumers),day_k] # log(density_current));
					end;
				end;
				
				**calculating acceptance probability;
				if log_likelihood_candidate >= log_likelihood_current then do;
				
					acceptance_probability = 1;
				end;
				else do;
				
					acceptance_probability = exp(log_likelihood_candidate - log_likelihood_current);
				end;
				
				**select candidate beta1 with probability equal to the acceptance probability, otherwise reject the candidate and keep the current beta1;
				acceptance_choice = j(1, 1, .);
				call randgen(acceptance_choice, "Uniform");
				beta_updated$1 = beta1_candidate*(acceptance_choice <= acceptance_probability) + beta1_current*(acceptance_choice > acceptance_probability);
			end;
			else do;
/* BW CHANGE 2025-10-22 */			
				**calculate the mean of the updated beta;
				*beta_mean = c2 * c1;
				beta_mean = matMultAB(c2, c1, alg);
				
				**calculate the standard deviation of the updated beta by taking the matrix square root of C2;
				*call eigen(c2_eigvals, c2_eigvecs, c2);
				*beta_std_dev = c2_eigvecs * diag(sqrt(c2_eigvals)) * c2_eigvecs`;
				tmp        = eigenDecomp(c2, alg);
				c2_eigvals = tmp$"eigenvalues";
				c2_eigvecs = tmp$"eigenvectors";
				free tmp;
				beta_std_dev = matMultADAt(c2_eigvecs, sqrt(c2_eigvals), alg);
				
				**calculate updated beta for this variable;
				**normal distribution with a mean of C2*C1 and a standard deviation of sqrt(C2);
				normals = j(num_coefficients_j, 1, .);
				call randgen(normals, "Normal");
				
				*beta_updated$var_j = beta_mean + beta_std_dev * normals;
				beta_updated$var_j = beta_mean + matMultAB(beta_std_dev, normals, alg);
			end;
		end;
		
		if symgetn("DEBUG") then print "end module: update_beta";
		return beta_updated;
	finish;

	**Updates conditional Ni for the first episodic food for each subject;
	**Negative for never-consumers and positive for consumers;
	**Uses complete conditional in Appendix A.2 of Bhadra, et al. (2020);
	start update_conni1(alpha1,
											consumer_probabilities,
											xbeta_u,
											never_consumer_covariate_matrix,
											episodic_indicator_matrices,
											recall_availability,
											num_subjects,
											num_recalls);
		
		if symgetn("DEBUG") then print "begin module: update_conni1";
		alg = symgetn("v143code") + 1;
		
		**probability that each subject reports no consumption days;
		prob_no_consumption_days = j(num_subjects, 1, 1);
		do day_k = 1 to num_recalls;
		
			xbeta_u_k = xbeta_u$day_k;
			
			prob_no_consumption_day_k = j(num_subjects, 1, 1);
			if any(recall_availability[,day_k]) then do;
			
				prob_no_consumption_day_k[loc(recall_availability[,day_k])] = 1 - cdf("Normal", xbeta_u_k[loc(recall_availability[,day_k]),1]);
			end;
			
			prob_no_consumption_days = prob_no_consumption_days # prob_no_consumption_day_k;
		end;
		
		**probability that each subject is a never-consumer given that they report no consumption;
		never_consumer_prob = 1 - consumer_probabilities;
		cond_never_consumer_prob = never_consumer_prob/(never_consumer_prob + consumer_probabilities # prob_no_consumption_days);
		
		**linear predictors for consumer probabilities;
/* BW CHANGE 2025-10-22 */
		*g_alpha = never_consumer_covariate_matrix * alpha1;
		g_alpha = matMultAB(never_consumer_covariate_matrix, alpha1, alg);

		**find subjects with with with no reported consumption;
		no_consumption_indicator = j(num_subjects, 1, 1);
		do day_k = 1 to num_recalls;
		
			episodic_indicator_matrix_k = episodic_indicator_matrices$day_k;
			no_consumption_indicator = no_consumption_indicator # (1 - episodic_indicator_matrix_k[,1]);
		end;
		
		**choose never-consumers based on conditional never-consumer probabilities;
		never_consumer_choice = j(num_subjects, 1, .);
		call randgen(never_consumer_choice, "Uniform");
		never_consumers = no_consumption_indicator & (never_consumer_choice <= cond_never_consumer_prob);
		
		**truncated normals;
		truncation_values = (2 # never_consumers - 1) # g_alpha;
		trunc_normal = truncated_normal(truncation_values);
		
		**update conni1;
		conni1_updated = j(num_subjects, 1, .);
		
		if any(^never_consumers) then do;
		
			conni1_updated[loc(^never_consumers)] = g_alpha[loc(^never_consumers)] + trunc_normal[loc(^never_consumers)];
		end; 
		
		if any(never_consumers) then do;
		
			conni1_updated[loc(never_consumers)] = g_alpha[loc(never_consumers)] - trunc_normal[loc(never_consumers)];
		end;
		
		if symgetn("DEBUG") then print "end module: update_conni1";
		return conni1_updated;
	finish;
	
	**Updates alpha1;
	**Uses complete conditional in Appendix A.3 of Bhadra, et al. (2020);
	start update_alpha1(alpha1_mean_prior,
											alpha1_covariance_prior,
											conni1,
											subject_weighting,
											never_consumer_covariate_matrix,
											num_never_consumer_coefficients);
		
		if symgetn("DEBUG") then print "begin module: update_alpha1";
		alg = symgetn("v143code") + 1;
		
		wt_covariate_matrix = (subject_weighting # never_consumer_covariate_matrix)`;
		
		**Calculate C1 (mean) in the complete conditional for alpha in Appendix A.3 of Bhadra, et al. (2020);
/* BW CHANGE 2025-10-22 */
		*c1_first_term = inv(alpha1_covariance_prior) * alpha1_mean_prior;
		*c1_second_term = wt_covariate_matrix * conni1;
		c1_first_term  = matMultAB(inv(alpha1_covariance_prior), alpha1_mean_prior, alg);
		c1_second_term = matMultAB(wt_covariate_matrix, conni1, alg);
		
		c1 = c1_first_term + c1_second_term;
		
		**Calculate C2 (variance) in the complete conditional for alpha in Appendix A.3 of Bhadra, et al. (2020);
		*wt_covariate_sq_sums = wt_covariate_matrix * never_consumer_covariate_matrix;
		wt_covariate_sq_sums = matMultAB(wt_covariate_matrix, never_consumer_covariate_matrix, alg);
		
		inverse_c2 = inv(alpha1_covariance_prior) + wt_covariate_sq_sums;
		c2 = inv(inverse_c2);
		
		**calculate the standard deviation of alpha (square root of C2);
/* BW CHANGE 2025-10-22 */	
		*call eigen(c2_eigvals, c2_eigvecs, c2);
		tmp        = eigenDecomp(c2, alg);
		c2_eigvals = tmp$"eigenvalues";
		c2_eigvecs = tmp$"eigenvectors";
		free tmp;
		
		*alpha1_std_dev = c2_eigvecs * diag(sqrt(c2_eigvals)) * c2_eigvecs`;
		alpha1_std_dev = matMultADAt(c2_eigvecs, sqrt(c2_eigvals), alg);
		
		**calculate updated alpha1;
		*alpha1_mean = c2 * c1;
		alpha1_mean = matMultAB(c2, c1, alg);
		
		normals = j(num_never_consumer_coefficients, 1, .);
		call randgen(normals, "Normal");
/* BW CHANGE 2025-10-22 */
		*alpha1_updated = alpha1_mean + alpha1_std_dev * normals;
		alpha1_updated = alpha1_mean + matMultAB(alpha1_std_dev, normals, alg);
		
		if symgetn("DEBUG") then print "end module: update_alpha1";
		return alpha1_updated;
	finish;
	
	**Updates consumer probabilities;
	start update_consumer_probabilities(never_consumer_covariate_matrix,
																			alpha1);
		
		if symgetn("DEBUG") then print "begin module: update_consumer_probabilities";
		alg = symgetn("v143code") + 1;
/* BW CHANGE 2025-10-22 */		
		*g_alpha = never_consumer_covariate_matrix * alpha1;
		g_alpha = matMultAB(never_consumer_covariate_matrix, alpha1, alg);
		
		consumer_probabilities = cdf("Normal", g_alpha);
		
		if symgetn("DEBUG") then print "end module: update_consumer_probabilities";
		return consumer_probabilities;
	finish;
	
	reset storage=WORK.mcmc_modules;

	store module=_all_;
quit;