Understanding sampling of parameters / local or model block variables

General
1

I just realized I do not understand how parameters are treated in Stan.

I built a simple physics model for a thrown object, and estimate the angle thrown based on time and distance. It is straightforward, and it works fine. (See code at end.)

The problem is that I don’t want to estimate the individual angles - for now, I only want to estimate the variance. But if I move the Calc_Angle variable to be a local variable in the model block, the sampler evidently isn’t exploring the space of possible values for those anymore, so the model can’t find a fit for the variances, and therefore doesn’t find a log probability between -2,2, and gives up.

2 questions:

  1. Am I understanding that correctly?

  2. Is there any way to tell the model to sample a parameter, but not include it in the model that is output?

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>[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.
  real<lower=0> g; // Gravity - 9.8 m/s on earth.
}

parameters {
  vector[N_i] Calc_Angle[N_p];
  real<lower=0> angle_tau;
  real<lower=0> overall_distance_tau;
}

model {
	angle_tau ~ gamma(4,1);
	overall_distance_tau ~ gamma(1,1);
  {  
  for (p in 1:N_p){
  Calc_Angle[p] ~ normal(angle_calc(measured_speed[p], measured_time[p], to_vector(rep_array(height_p[p],N_i)),g), angle_tau);
  measured_distance[p] ~ normal(dist_calc(measured_speed[p], Calc_Angle[p], to_vector(rep_array(height_p[p],N_i)),g), overall_distance_tau);
    }
  }
}
2

The Stan program you write is only describing the calculation of the density and its gradient with respect to the parameters. That is it. Then you pick a point in parameter space (initial values) and HMC explores from that point by moving through parameter space. Even though calc_angle appears in your model on the LHS of an expression, you are never using Stan code to calculate its value, so this line:

Calc_Angle[p] ~ normal(angle_calc(measured_speed[p], measured_time[p], to_vector(rep_array(height_p[p],N_i)),g), angle_tau);

Is really syntactic sugar (that I hate) for:

target += normal_lpdf(Calc_Angle[p] | angle_calc(measured_speed[p], measured_time[p], to_vector(rep_array(height_p[p],N_i)),g), angle_tau);

If you don’t add some way of calculating Calc_angle then it starts out as either (I don’t remember which) an arbitrary value or NaN or 0. So the log density calculation is never completed successfully and initialization is rejected. In any case, not what you want.

1 Like
3

In RStan, yes. You can call stan with pars = "Calc_Angle", include = FALSE. The draws won’t be included in the output, although you can still get the posterior means. We need to do this for Stan generally.

1 Like