How to assign Different Contraints for Variables sampled from a Multivariate Normal?

Modeling
,
1

Hey all,

I’m trying to fitting a model with three free parameters from a MVN (full model code at the end, parts to emphasize the problem). The Problem is, that one Parameter is on the logit scale, whereas others are on are normal. Thus, I have to assign different constraints for the hyper and the subject parameters for the parameters, if I want to them to sample from a MVN, with reasonable results).

I managed the hyper par contraints with the following code:

data 
{ 
vector[J] L;  // vector of three lower boundaries
}
parameters
{    
vector<lower=0>[J] hyper_pars_raw;
}
model
{
 hyper_pars = L + hyper_pars_raw;
}

This seems to work fine, but I have problems to apply the same Solution with the subject parameters. I did it this way:

data {

vector[J] L;   // vector of three lower boundaries

}

parameters {
 vector<lower=0>[J] subj_pars_raw[N];
}

// Apply Bounds for Subject Pars

transformed parameters
 {
      for (k in 1:N) // is the Number of Subjects.
          {
              subj_pars[k,1] = L[1] + subj_pars_raw[k,1];
              subj_pars[k,2] = L[2] + subj_pars_raw[k,2];
              subj_pars[k,3] = L[3] + subj_pars_raw[k,3];
           }
}

model
{

 target += sum(hyper_pars);

// priors for hyper parameters

 Omega ~ lkj_corr(1);
  sigma ~ gamma(1,0.01);

hyper_pars[1] ~ normal(20,5);
hyper_pars[2] ~ normal(20,5);
hyper_pars[3] ~ normal(0,10);

// Sample Parameters

subj_pars[,1:3] ~ multi_normal(hyper_pars, Sigma);

  for(i in 1:N)
   {
     // draw data from probabilities determined by parms
      count[i,]  ~ multinomial(probs[i,]);*
   }

}

I got now the following error, which I can’t relate:


Chain 1: Rejecting initial value:
Chain 1:  Error evaluating the log probability at the initial value.
Chain 1: Exception: multinomial_lpmf: Probabilities parameter is not a valid simplex. 
sum(Probabilities parameter) = nan, but should be 1
(in 'model567430963fc4_ComplexSpan_LKJ' at line 114)

I also tried to define two varibales as constrained vector and one as unconstrained vector, to bound them together with append_col. But as the MVN needs a vector and can’t handle Matrices, this also didn’t worked.
I transformed it with to_vector, but got a dimension error then.

Here is the full model code:

data {
  int <lower=0> N;  // number of subjects
  int <lower=0> K;  // categories 
  int <lower=0> J;  // Dims of Cov Matrix
  vector[J] L;
  int R[K];         // number of responses per category
  int count[N,K];   // observed data
  real scale_b;     // set scaling for background noise
}

parameters {
  // subject parameters
  
  corr_matrix[J] Omega;
  vector<lower=0>[J] sigma;
  vector<lower=0>[J] hyper_pars_raw;
  vector<lower=0>[J] subj_pars_raw[N];
  
  
}

transformed parameters{

  vector[J] subj_pars[N]; 
  vector[J] hyper_pars;
  cov_matrix[J] Sigma; 
  
  // Transform f Parameter
  
  real f[N] = inv_logit(subj_pars[,3]);
  real mu_f = inv_logit(hyper_pars[3]);
  
  // activations
  vector[K] acts[N];
  real SummedActs[N];
  
  // probabilities
  
  vector[K] probs[N];
  
  
  // Transformations
  
   hyper_pars = L + hyper_pars_raw;

  
  for (k in 1:N)

{
  subj_pars[k,1] = L[1] + subj_pars_raw[k,1];
  subj_pars[k,2] = L[2] + subj_pars_raw[k,2];
  subj_pars[k,3] = L[3] + subj_pars_raw[k,3];


}



  
  // Transform Sigma
  Sigma = quad_form_diag(Omega, sigma); 
  
  // loop over subjects to compute activations and probabilites
  
  for (i in 1:N){
    acts[i,1] = scale_b + subj_pars[i,1] + subj_pars[i,2]; // Item in Position
    acts[i,2] = scale_b + subj_pars[i,1];        // Item in Other Position
    acts[i,3] = scale_b + f[i]* (subj_pars[i,1]+subj_pars[i,2]);// Distractor in Position
    acts[i,4] = scale_b + f[i]*subj_pars[i,1]; // Distractor in other Position
    acts[i,5] = scale_b; // non presented Lure
    
    SummedActs[i] = R[1] * acts[i,1] + R[2] * acts[i,2] + R[3] * acts[i,3]+ R[4] * acts[i,4]+ R[5] * acts[i,5];
    
    probs[i,1] = (R[1] * acts[i,1]) ./ (SummedActs[i]);  
    probs[i,2] = (R[2] * acts[i,2]) ./ (SummedActs[i]);
    probs[i,3] = (R[3] * acts[i,3]) ./ (SummedActs[i]);
    probs[i,4] = (R[4] * acts[i,4]) ./ (SummedActs[i]);
    probs[i,5] = (R[5] * acts[i,5]) ./ (SummedActs[i]);
    
    
  }
  
}


model{
  
  target += sum(hyper_pars);

  // priors for hyper parameters
  
  Omega ~ lkj_corr(1);
  sigma ~ gamma(1,0.01);
  
  hyper_pars[1] ~ normal(20,5);
  hyper_pars[2] ~ normal(20,5);
  hyper_pars[3] ~ normal(0,10);
  
  // Sample Parameters
 
  subj_pars[,1:3] ~ multi_normal(hyper_pars, Sigma);

  // Loop over subjects
  for(i in 1:N){
    
    // draw data from probabilities determined by MMM parms
    count[i,]  ~ multinomial(probs[i,]);
  }
}

Thanks for your time! Hopefully someone has similiar issues!

1 Like
2

@jangoe indicated in a message to staff the problem was solved. We are happy for you and it woudl be great, if you could at least outline what the issue was so that others can learn :-)

Best of luck with your project!

2 Likes
3

I solved the issue by letting all parameters unconstrained - to avoid the underflow I defined a functions which randomly samples reasonable initial values from a uniform distribution for the subject parameters. The model runs fast and stable with this solution!

1 Like