Median of a vector

Has anyone implemented a function to calculate the median of a vector? I’ve given it a try myself, but ran into lots of complications with types (indices needing to be converted to floats to divide them by 2 and then back to integers to serve as indices again). So I wanted to check if perhaps there’s already a working implementation out there.

Thanks so much!


As of Stan 2.27, we have a quantile function: 6.6 Reductions | Stan Functions Reference

Which you would use as:

real median = quantile(your_vector, 0.5);

That’s exactly what I was looking for, thank you!


1 Like

Actually, I ran into some problems using this function. I need to take the median of some parameters, not of some data. So I got the error:

Ill-typed arguments supplied to function 'quantile':
(vector, real)
Available signatures:
(data vector, data real) => real
  The first argument must be data-only. (Local variables are assumed to
  depend on parameters; same goes for function inputs unless they are marked
  with the keyword 'data'.)
(data vector, data array[] real) => array[] real
  The first argument must be data-only. (Local variables are assumed to
  depend on parameters; same goes for function inputs unless they are marked
  with the keyword 'data'.)
(data row_vector, data real) => real
  The first argument must be row_vector but got vector
(data row_vector, data array[] real) => array[] real
  The first argument must be row_vector but got vector
(data array[] real, data real) => real
  The first argument must be array[] real but got vector
(Additional signatures omitted)

Is there a way I can use the function on a vector of parameters?

The problem is that the first argument is declared as data, which means it only works when applied to expressions that only involve data or transformed data variables. Presumably it’s not implemented for parameter vectors because it’s not differentiable everywhere.

Sample medians are tricky if there is an even number of elements in the vector. So you could write a median function as follows.

bool is_even(int n) {
  return 2 * (n / 2) == n;
real median(vector x) {
  int N = cols(y);
  if (N == 0) reject("median(x): x must not be empty");
  if (N == 1) return x[1];
  vector y = sort_asc(x);
  int mid = N / 2;
  if (is_even(N)) {
   return (y[mid] + y[mid + 1]) / 2;
  } else {
    return y[mid + 1];

I haven’t tested, but something like that should work. Just don’t expect to be able to do anything with median(x) in the model. The derivative you’ll get is just w.r.t. the returned value, not the sort.

If your program allows moving the median calculation down to the generated quantities block, that’s a better solution.

Thanks so much, Bob. I was planning to calculate the median of a set of parameters in the transformed parameters section, basically so that I can center them by their median for identification purposes.

If I then use those centered parameters in the model block, do you think I will run into issues?

The context is that this is an ordinal Item Response Theory model, and it seemed quite intuitive to me to identify the cutpoints by centering them at their median. So that the mean of all the “center” cutpoints across all questions is forced to be zero.

Not sure if it makes a difference, but these parameters are already in ordered vectors. So I actually don’t need the sort step, and I could use the rest of your function to return what I do need: the middle value if length is uneven, and the mean of the two middles if length is even.

In that case, can I use the output of this median function, applied to parameters, to transform my parameters and then use those transformed versions in the model block?

I think it makes all the difference that the parameters are in ordered vectors. You don’t need @Bob_Carpenter’s nifty function, you just need to pull out the middle element or the mean of the middle 2 elements.

Alternatively, for an odd number of elements you can enforce a median-zero ordered vector by concatenating a negative ordered vector, zero, and a positive ordered vector. For an even number of elements, you can concatenate a negative ordered vector with maximum element m, -m, and -m + v where v is a positive ordered vector.

1 Like

Thanks, Jacob! I thought I’d use Bob’s function anyway, because it does exactly that if we skip the sort statement. Since I have many of these vectors and I won’t know how many elements they have, I do need a function like this to get at the median of each of them.

To clarify: I don’t want the median of these to become zero, I want to find the median so that I can use the average of all those medians in a centering step. Like so:

  //calculate average median threshold across items
  for(j in 1:J){
      mean_beta_mid += median( segment(beta_uncentered[j], 1, Kmin1[j]));
  mean_beta_mid = mean_beta_mid / J;
  //center betas on this average
  for(j in 1:J){
    beta[j] = beta_uncentered[j] - mean_beta_mid;

So my remaining issue is whether or not this will break the derivative. When I compile this in pedantic mode, I do get a warning about a control flow statement depending on a parameter.

1 Like

Makes sense! Since you want to center by an average of medians rather than by a single median the concatenation tricks I mentioned won’t work, as you’ve rightly pointed out. However, there might be some models where it could be sensible to parameterize your models in terms of these zero-median vectors plus the true median. I cannot say whether or not this would actually ever be a good idea.

If it matters to you, because the vectors are ordered you can find the median more efficiently than in the more general implementation that @Bob_Carpenter provided. In particular, you just have to check whether the length is odd or even, and then return either the N/2 + 0.5 th element (write a function to return the smallest integer greater than or equal to N/2; Stan’s strong typing of integers will be a bit annoying here, but it can be dealt with) or the mean of the N/2 and N/2 + 1 elements.

Thanks! To report back briefly on the downstream effects of using these median-centered parameters: it hurt convergence a lot. I’m seeing Rhats of 2 and over on key model parameters.

I’m switching back to centering on the mean, instead. I got good convergence with that approach, so it’s definitely the median-centering step that ruined convergence.