How to set hard constraints on a transformed parameter within a hierarchical model?

Hi there, I have developed a hierarchical model of behavioural data in a psychophysical experiment. The model defines a grand average prior on parameters ‘bias’, ‘sens’, and ‘n_exp’. The model also defines per-session and per-subject variations to these parameters, which is where the hierarchical aspect of the model comes in. See code below.

So for example, the N_EXP parameter for subject 1 on session 1 is set as:
N_EXP = n_exp + <deviation for subject 1> + <deviation for session 1>
In the code below these deviations are held in the matrices b_sess and b_subj

One problem I’m having is that the variation in the parameters across sessions and subjects occasionally sets this N_EXP value to be less than zero, even though the global parameter n_exp was constraints to be positive. This is a problem when computing the likelihood because the contrastLeft or contrastRight values (which are raised to the power N_EXP) are sometimes 0. This makes the term 0^(negative number) which is undefined.

Therefore is there a way to constrain N_EXP transformed parameter like this, given that it is computed differently subject to subject and session to session?

Thank you

data {
	int<lower=1> numTrials;
	int<lower=1> numSessions;
	int<lower=1> numSubjects;
	int<lower=1,upper=numSessions> sessionID[numTrials];
	int<lower=1,upper=numSubjects> subjectID[numTrials];
	vector<lower=0,upper=1>[numTrials] contrastLeft;
	vector<lower=0,upper=1>[numTrials] contrastRight;
	int<lower=1,upper=3> choice[numTrials]; // 1=Left, 2=Right, 3=NoGo
	int<lower=1,upper=numSubjects> subjID_session[numSessions]; //the subject ID for each session in the dataset
parameters {
	//global parameters
	vector[2] bias; 
	vector<lower=0>[2] sens;
	real<lower=0> n_exp;
	//per session deviations 
	vector<lower=0>[5] sd_sess;
	matrix[5,numSessions] z_sess; //standard normal variable used to draw from the covariance matrix
	cholesky_factor_corr[5] rho_sess; //correlations of deviations
	//per subject deviations
	vector<lower=0>[5] sd_subj;
	matrix[5,numSubjects] z_subj; 
	cholesky_factor_corr[5] rho_subj; 
transformed parameters {
	vector[3] logOdds[numTrials]; // ln pL/pNG, ln pR/pNG, ln pNG/pNG
	matrix[5,numSessions] b_sess;
	matrix[5,numSubjects] b_subj;
	//draw samples of sess and subj deviations, according to the covariance structure in rho_ & sd_
	b_sess = diag_pre_multiply(sd_sess, rho_sess) * z_sess;
	b_subj = diag_pre_multiply(sd_subj, rho_subj) * z_subj;

		//temp variables
		real BL;
		real BR;
		real SL;
		real SR;
		real N_EXP;
		//compute (non)linear model
		for (n in 1:numTrials)
			BL = bias[1] + b_sess[1,sessionID[n]] + b_subj[1,subjectID[n]];
			BR = bias[2] + b_sess[2,sessionID[n]] + b_subj[2,subjectID[n]];
			SL = sens[1] + b_sess[3,sessionID[n]] + b_subj[3,subjectID[n]];
			SR = sens[2] + b_sess[4,sessionID[n]] + b_subj[4,subjectID[n]];
			N_EXP= n_exp + b_sess[5,sessionID[n]] + b_subj[5,subjectID[n]];

			logOdds[n][1] = BL + SL*contrastLeft[n]^N_EXP;
			logOdds[n][2] = BR + SR*contrastRight[n]^N_EXP;
			logOdds[n][3] = 0;
model {
	//priors on global parameters
	bias ~ normal(0, 2);
	sens ~ normal(5, 2);
	n_exp ~ normal(0.5,0.25);
	//make z_std_normal be standard normal (non centred parameterisation)
	to_vector(z_sess) ~ normal(0, 1);	
	to_vector(z_subj) ~ normal(0, 1);	
	//prior on the variation of the per-session deviations
	sd_sess ~ cauchy(0,1);
	sd_subj ~ cauchy(0,1);
	//prior on the cholesky factor of the covariance matrix
	rho_sess ~ lkj_corr_cholesky(2.0); //penalises extreme correlations between the deviations
	rho_subj ~ lkj_corr_cholesky(2.0); //penalises extreme correlations between the deviations

	for (n in 1:numTrials) {
		choice[n] ~ categorical_logit( logOdds[n] );
generated quantities {
	corr_matrix[5] corr_sess;
	corr_matrix[5] corr_subj;
	vector[numTrials] log_lik;
	//write correlation matrix
	corr_sess = rho_sess * rho_sess';
	corr_subj = rho_subj * rho_subj';

	//write loglik
	for (n in 1:numTrials){
		log_lik[n] = categorical_logit_lpmf(choice[n] | logOdds[n] );

I realize this is an old post, but I have the same question.

Does anyone have an answer to this question?

Thank you,