Problem when converting to a ragged array model (unsure why - syntax? typo?)

Modeling
1

I’m attempting to convert a model that uses rectangular data into one that can deal with ragged arrays.

Step 1 is to rewrite it using non-ragged data entered in the new format and just replicate the results, based on the standard trick for looping through an arrays indexed by individual - but Stan is telling me that my code is wrong, I suspect, (Error evaluating the log probability at the initial value.) and I’ve failed at seeing where the difference is.

I am inputting the same data, except that I’m not transforming it into a matrix in the ragged case - but I can’t figure out how to look at any diagnostics, and stan won’t display anything I print() in the parameters and model blocks because it terminates before sampling. :(

I tried replicating with a simpler model, but in that case it works, so I’m unsure what is different. Am I being silly, and missing something obvious, or should these two models not return the same result given the same data?

Non-ragged: (works)

functions {
vector dist_calc(vector speed, vector angle, vector start_height, real g){
  vector[size(to_array_1d(speed))] v_h;
  v_h = speed .* sin(angle*pi()/180);
  return speed .* cos((angle)*pi()/180) .* ((v_h + sqrt(v_h .* v_h + 2*g*start_height)) / (g));
  }
vector angle_calc(vector speed, vector time, vector start_height, real g){
  // (180/pi)* asin( ((1/2) * g * measured_time^2 - height_p[Thrower_i])/(measured_speed * measured_time))
  return asin(((g/2*(time .* time)-start_height)) ./ (speed .* time))*(180/pi());
  } 
}

data {
  //Data Element Parameters;
  int<lower=1> N_i; //Number of input cases per person
  int<lower=1> N_p; //Number of throwers
  int<lower=1> Thrower_i[N_i*N_p]; //Which thrower threw case i
  vector<lower=0,upper=90>[N_i] throw_angle_intent[N_p]; //Angle Intended for case i - but varies based on the error in each throw
  vector<lower=0>[N_i] measured_speed[N_p]; //Speed Measured for case i in M/s
  vector<lower=0>[N_i] measured_time[N_p]; //Time Measured for case i from throw until landing
  vector<lower=0>[N_i] measured_distance[N_p]; //Distance Measured for case i in meters
  vector<lower=0>[N_p] height_p; // Person P's height.

  // Model Parameter elements:
  real<lower=0> g; // Gravity - 9.8 m/s on earth.
}


transformed data {
}

parameters {
  vector[N_i] Angle_Variance[N_p];
  real<lower=0> angle_tau;
  vector[N_p] personal_angle_bias;
  real<lower=0> dist_calc_tau;
  } //End Parameters

transformed parameters {
}

model {
	angle_tau ~ gamma(4,1);
	dist_calc_tau ~ gamma(10,1);
  {  
  for (p in 1:N_p){
  vector[N_i] Calc_Angle;
  Calc_Angle = angle_calc(measured_speed[p], measured_time[p], to_vector(rep_array(height_p[p],N_i)),g); //This is determined by kinematics.
  Angle_Variance[p] ~ normal(0,5);
  throw_angle_intent[p] ~ normal(Angle_Variance[p] + Calc_Angle,0.001);
  personal_angle_bias[p] ~ normal(mean(Angle_Variance[p]),0.001);
  measured_distance[p] ~ normal(dist_calc(measured_speed[p], Calc_Angle, to_vector(rep_array(height_p[p],N_i)),g), dist_calc_tau);
  }
  }
}

Ragged version: (fails)

functions {

vector dist_calc(vector speed, vector angle, vector start_height, real g){
  vector[size(to_array_1d(speed))] v_h;
  v_h = speed .* sin(angle*pi()/180);
  return speed .* cos((angle)*pi()/180) .* ((v_h + sqrt(v_h .* v_h + 2*g*start_height)) / (g));
  }

vector angle_calc(vector speed, vector time, vector start_height, real g){
  // (180/pi)* asin( ((1/2) * g * measured_time^2 - height_p[Thrower_i])/(measured_speed * measured_time))
  return asin(((g/2*(time .* time)-start_height)) ./ (speed .* time))*(180/pi());
  } 
}
data {
  //Data Element Parameters;
  int<lower=1> N_i_t; //Number of total throws
  int<lower=1> N_p; //Number of throwers
  int<lower=1> Thrower_i[N_i_t]; //Which thrower threw case i
  vector<lower=0,upper=90>[N_i_t] throw_angle_intent; //Angle Intended for case i - but varies based on the error in each throw
  vector<lower=0>[N_i_t] measured_speed; //Speed Measured for case i in M/s
  vector<lower=0>[N_i_t] measured_time; //Time Measured for case i from throw until landing
  vector<lower=0>[N_i_t] measured_distance; //Distance Measured for case i in meters
  vector<lower=0>[N_p] height_p; // Person P's height.
  // Model Parameter elements:
  real<lower=0> g; // Gravity - 9.8 m/s on earth.
}


transformed data {
int throws[N_p];
throws = rep_array(0, N_p);
for (i in 1:N_i_t){
	throws[Thrower_i[i]] = throws[Thrower_i[i]]+1; //How many throws each player had.
}
print(throws);
print(measured_speed);
print(measured_time);
print(height_p);
}

parameters {
  vector[N_i_t] Angle_Variance;
  real<lower=0> angle_tau;
  vector[N_p] personal_angle_bias;
  real<lower=0> dist_calc_tau;
  } //End Parameters

transformed parameters {
}

model {
	angle_tau ~ gamma(4,1);
	dist_calc_tau ~ gamma(10,1);

  {
  int position;
  position = 1;
  
  for (p in 1:N_p){
  vector[throws[p]] measured_speed_p;
  vector[throws[p]] measured_time_p;
  vector[throws[p]] measured_distance_p;
  vector[throws[p]] Angle_Variance_p;
  vector[throws[p]] throw_angle_intent_p;

  vector[throws[p]] Calc_Angle;
  
  measured_speed_p = to_vector(segment(measured_speed,position, throws[p]));
  measured_time_p = to_vector(segment(measured_time,position, throws[p]));
  measured_distance_p = to_vector(segment(measured_distance,position, throws[p]));
  Angle_Variance_p = to_vector(segment(Angle_Variance,position, throws[p]));
  throw_angle_intent_p = to_vector(segment(throw_angle_intent_p,position, throws[p]));
  
  Calc_Angle = angle_calc(measured_speed_p, measured_time_p, to_vector(rep_array(height_p[p],throws[p])),g); //This is determined by kinematics.
  Angle_Variance_p ~ normal(0,angle_tau);
  throw_angle_intent_p ~ normal(Angle_Variance_p + Calc_Angle,0.001);
  personal_angle_bias[p] ~ normal(mean(Angle_Variance_p),0.001);
  measured_distance[p] ~ normal(dist_calc(measured_speed_p, Calc_Angle, to_vector(rep_array(height_p[p],throws[p])),g), dist_calc_tau);
  
  position = position +  throws[p]; //Move the offset to be correct for the next player.
  }
  }
}
2

Sorry about that. The current RStan (and maybe other interfaces) are swallowing warning messages during initialization (and prints during rejections). This is motivating us to get 2.16 sooner than we would otherwise, as nobody can debug their complicated models any more.

You really don’t want to redefine those data constants as local variables as it will create a bunch of autodiff variables.

We use scale, not precision, so that scale of 0.001 for a normal is probably an error. If you’re trying to just lock things down, it’s hard for Stan to adapt that small, but you can help by rescaling manually.

You can use compound declare define andyou can also use the slicing notation, which should all be easier to read than to_vector(segiment(...)).

int start = 1;
for (p in 1:N_p) {
  int num_throws = throws[p];
  int end = start + num_throws;
  ... then use, e.g, measured_time[start:end] instead of defining new variabels
 ...
}
1 Like
3

I have followed your suggestion about removing the local variables, and at least it compiles more quickly - though I have no idea if it works, since I get the same error, and I tried redoing things to not need a variance that is effectively 0, which should also help - but it works fine in the non-ragged model, so…

I’m now trying to run this in cmdstan instead, which presumably will display the error messages - but running into problems there. (See my other post.)

4

Hi Bob,
Now that I am running cmdstan, I’m still not seeing anything I can use to figure out what is wrong - no additional error messages are showing up.

I also still cannot seem to print anything from the transformed variables block :(

Any suggestions for what I should be looking for?

5

@davidmanheim Do you have the updated ragged model with the segments changed to slices and the normal(x, 0.001) constraints removed? I don’t mind looking over it. The error wasn’t obvious to me in the last model, but maybe if the pieces are shuffled around a little it’ll be easier.

1 Like
6

Thanks - but I figured it out, eventually.

This part:

  throw_angle_intent_p ~ normal(Angle_Variance_p + Calc_Angle,0.001);
  personal_angle_bias[p] ~ normal(mean(Angle_Variance_p),0.001);

was giving the model a ton of trouble, because it was sampling for each independently, but in fact the second should be treated as a deterministic outcome of the first.

Also, for Bob’s suggested fix, I needed to be careful - the indexing given had an off-by-1 error that stan didn’t notice, but would have broken things anyways, if I hadn’t noticed it while looking for other issues.
I needed to fix the offset: end = start + num_throws -1.

1 Like
7

Thanks for following up. Looks like the fix won’t be in until 2.16, which is in the works.