# Returning quantiles (or weighted quantiles) in generated quantities block

**URL:** <https://discourse.mc-stan.org/t/returning-quantiles-or-weighted-quantiles-in-generated-quantities-block/28787>\
**Category:** Modeling\
**Tags:** rstan\
**Created:** [September 8, 2022, 6:48pm UTC](https://discourse.mc-stan.org/t/returning-quantiles-or-weighted-quantiles-in-generated-quantities-block/28787 "2022-09-08T18:48:45Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![AFH](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/afh/32/15068_2.png) [@AFH](https://discourse.mc-stan.org/u/AFH)\
**Post date:** [September 8, 2022, 6:48pm UTC](https://discourse.mc-stan.org/t/returning-quantiles-or-weighted-quantiles-in-generated-quantities-block/28787/1 "2022-09-08T18:48:45Z")

</div>

Hi, I’m running a complex Stan model via RStan. It estimates population densities over space and time. I want to calculate some summary statistics of these populations (like center of spatial distribution, edge of spatial distribution, etc) within Stan as generated quantities, so they have their own posterior distributions, rather than calculating them post hoc in R.

Some of these are calculated as weighted quantiles. For example, the center of spatial distribution is the median of observations across space, weighted by the size (abundance) of each observation. The edge is usually calculated a 0.05 or 0.95 quantile of the same.

I know a quantile function has recently been implemented in Stan, but I’m having trouble getting it to work—the [function reference link is broken](https://mc-stan.org/docs/functions-reference/phi-function.html) and I haven’t totally parsed the long GH issues discussing its development.

I tried writing a quick model to explore how the quantile function works, but can’t even get that to work. So, this is a two-part question:

1. HOW DOES QUANTILE() WORK IN STAN?

My toy model:

```nohighlight
data {
real probs;
int np;
int ny;
}

parameters {
  matrix<lower=0>[np, ny] N;  
}

model{
  for(p in 1:np){
    for(y in 1:ny){
      N[p,y] ~ normal(10000, 3000);
    }
  }
}

generated quantities {
    real quant_out[ny]; 
    for(y in 1:ny){
quant_out = quantile(N[,y], probs); 
    }
}

```

Returns this error that seems to imply that quantile() doesn’t return anything:

```nohighlight
A returning function was expected but an undeclared identifier 'quantile' was supplied.

```

And for extra credit…  
2. HAS ANYONE WRITTEN A USER-DEFINED FUNCTION FOR CALCULATING WEIGHTED QUANTILES IN STAN?

This function would take as arguments a real probability / quantile, a vector or 1D array of data, and a vector or 1D array of weights (same dimensions as the data). I understand there are multiple ways to calculate weighted quantiles, and I’m not too picky about the exact method, since I’ll explore sensitivity to it regardless.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [September 20, 2022, 10:43pm UTC](https://discourse.mc-stan.org/t/returning-quantiles-or-weighted-quantiles-in-generated-quantities-block/28787/2 "2022-09-20T22:43:37Z")

</div>

here’s the link to the function reference for the standard normal cdf `Phi()` if that’s what you wanted:

> **[3.12 Probability-related functions | Stan Functions Reference](https://mc-stan.org/docs/functions-reference/Phi-function.html)**
>
> Reference for the functions defined in the Stan math library and available in the Stan programming language.

We also have `_cdf` and `_lcdf` for cdfs and log cdfs.

> [@AFH](#):
>
> 1. HOW DOES QUANTILE() WORK IN STAN?

I was about to say there wasn’t a `quantile()` function in Stan, but apparently someone added one; the doc is here:

- [6.6 Reductions | Stan Functions Reference](https://mc-stan.org/docs/functions-reference/reductions.html)

These only work for data, not for parameters, so your program isn’t going to compile. You need at least Stan 2.27, which means you are going to have to use `cmdstan`, `cmdstanr`, or `cmdstanpy` (or maybe `PyStan`, but `rstan` is out of date on CRAN and is only up to 2.26 on GitHub (this is due to CRAN’s restrictive policies on package size and dependencies and the Stan project not organizing itself around CRAN releases).

And we also don’t have inverse CDF functions implemented other than for the normal distribution. But given that N[p, y] is a truncated `normal(1e4, 3e3)`, you can compute the actual quantile for the density, not the empirical one from your sample.

Also, if you can get parameters closer to unit scale, Stan will work a lot better. We don’t have a way to supply both a lower bound and an offset and multiplier. So you’re probably better off implementing this as:

```no-highlight
parameters {
  matrix<lower=0>[np, ny] N_std;
}
transformed parameters {
  matrix<lower=0>[np, ny] N = 3e3 + 1e4 * N_std;
}
model {
  to_vector(N_std) ~ std_normal();
}

```

Stan samples the parameters, so this puts the sampler back on the centered unit scale.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [September 20, 2022, 10:44pm UTC](https://discourse.mc-stan.org/t/returning-quantiles-or-weighted-quantiles-in-generated-quantities-block/28787/3 "2022-09-20T22:44:57Z")

</div>

> [@AFH](#):
>
> 1. HAS ANYONE WRITTEN A USER-DEFINED FUNCTION FOR CALCULATING WEIGHTED QUANTILES IN STAN?

I haven’t seen one, but then I don’t know what a weighted quantile is. If you can write down the definition, it shouldn’t be too hard to define. But the problem is that it’s probably not going to do what you want, which is why we only have quantiles implemented for data. The problem with selecting and branching is that we can’t autodiff through those choices, so you wind up disconnecting the densities.

---

<div class="post-metadata">

**Author:** ![AFH](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/afh/32/15068_2.png) [@AFH](https://discourse.mc-stan.org/u/AFH)\
**Post date:** [September 21, 2022, 12:20am UTC](https://discourse.mc-stan.org/t/returning-quantiles-or-weighted-quantiles-in-generated-quantities-block/28787/4 "2022-09-21T00:20:59Z")

</div>

Thanks so much for all of this valuable explanation, I really appreciate it. I got approximately the function I wanted since posting this (apologies that there are some use-case specific things in here about indexing over patches and years, as I’m modeling changing species distributions). The note that `quantile()` only works for data explains a lot! Your reply gives me yet another reason to transition this implementation from `rstan` to `cmdstanr`.

```nohighlight

  // cumulative sum function
  // slightly different than the built-in Stan one -- this takes a 1D array of numbers and returns a single real number that sums elements 1:x of that array
  real csum(real[] to_sum, int x){
    return(sum(to_sum[1:x])); 
  }
  
  // function to calculate range quantiles
  
  real calculate_range_quantile(int np, vector patches, real[] dens_by_patch, real quantile_out){
    vector[np] csum_dens;
    vector[np] quant_p; 
    real diff_p[np]; 
    real quant_position; 
    real dens_p[np];
    
    for(i in 1:np){
      csum_dens[i] = csum(dens_by_patch[1:np], i); // calculate cumulative sum of density along each patch 
    }
    
    quant_p = csum_dens / sum(dens_by_patch[1:np]); // turn in to proportions, i.e., quantiles 
    
    for(i in 1:np){
      diff_p[i] = fabs(quant_p[i] - quantile_out); // minimizes absolute distance between estimated and desired quantile 
    }
    
    // DOESN'T DEAL WITH TIES! (but they seem unlikely with only 10 patches)
    
    // probably a more elegant way to do this... 
    for(i in 1:np){
      if(diff_p[i] == min(diff_p)) {
        quant_position = patches[i];
      }
    }
    return(quant_position); 
    
  } // close quantile function

```