Incomplete beta as prior

Modeling
1

Hi,

In the mean-shape parameterization of the beta distribution we have two parameters mu and theta. When theta = 2 every probability from zero to one is equally like. When theta < 2, the distribution is so dispersed that extreme probabilities near 0 and 1 are more likely than the mean.
Therefore, I want to make sure that theta >= 2.

In my model I introduce inv_theta in the “transformed parameters”-block. Logically, I need to make sure that inv_theta <= 0.5.

Simplified model:

...

parameters { 
  ...
  real<lower=0, upper=0.5> inv_theta;
}

transformed parameters {
  real<lower=2> theta  = 1.0 / inv_theta; 
}

model { 
  // prior
  inv_theta ~ exponential(1); 
  ....
}

Is this correctly done? Or should I make sure that inv_theta is drawn from a prior that generates only values between 0 and 0.5?
To my understanding there is no such distribution except the incomplete beta-distribution like incomplete_beta(0.5, a, b). But it’s not possible to use the incomplete beta-distribution as a prior, correct?

2

Assuming that a parameter is declared in the parameters block, appears on the left hand side of exactly one prior statement, and is not used to compute transformed parameters with prior statements of their own, then the realized prior will correspond exactly to the prior statement truncated by any constraints in the parameter declaration.

3

[edit: wasn’t done writing]

Note quite, because the following statement isn’t quite right:

If you want to avoid poles at 0 or 1, you have to make sure that

\min(\mu, 1 - \mu) \cdot \theta \geq 1.

or, equivalently

\theta \geq 1 / \min(\mu, 1 - \mu).

In Stan, that’s

parameters {
  real<lower=0, upper=1> mu;
  real<lower=1 / min(mu, 1 - mu)> theta;
}

If you want to work in terms of \phi = 1\theta, then we know that

1 / \phi \geq 1 / \min(\mu, 1 - \mu)

or equivalently,

\min(\mu, 1 - \mu) > \phi, or in Stan.

parameters {
  real<lower=0, upper=1> mu;
  real<upper = min(mu, 1 - mu)> inv_theta;
}
transformed parameters {
  real theta = inv(inv_theta);
}

Then you can put whatever kind of prior you want on inv_theta and it will be truncated appropriately to avoid poles at 0 or 1.

If you have a prior p(\phi), then it gets truncated from the declaration to

\displaystyle p(\phi) = \dfrac{\text{exponential}(\phi \mid 1)}{\int_0^{\min(\mu, 1 - \mu)} \text{exponential}(x \mid 1) \, \text{d}x}.

The way truncated declarations work in Stan, this happens automatically because the 1 in exponential(1) is a constant. If it wasn’t constant, you’d have to put the truncation in the distribution statement for inv_theta.

We had no idea what to call this and called it mean/precision. Technically, \theta is the sum of two shape parameters (beta terminology) or concentration parameters (Dirichlet terminology).

3 Likes
4

Thank you for the detailed answer. It’s very informative.

I have a final question that focuses more on the modelling itself.
In my model mu, inv_theta and theta are not reals, but they are defined as a matrix.

In this setting the constaint would look as follows:

min(\mu_{n,m}, 1− \mu_{n,m}) > \phi_{n,m}

If we write the below, does Stan works as intended?

parameters {
   matrix<lower=0, upper=1> mu;
   matrix<upper = min(mu, 1 - mu)> inv_theta;
}
transformed parameters {
   matrix theta = inv(inv_theta);
}

Thanks

5

Not quite. I think we can only accept arrays as lower and upper bounds. And you’ll need bounds. But what you can do is this:

data {
  int<lower=0> M, N;
  ...
parameters {
  vector[M * N] mu_vec;
  vector<upper = min(mu, 1 - mu)>[M * N] inv_theta_vec;

model {
   matrix[M, N] theta = inv(to_matrix(inv_theta_vec, M, N));
6

Thank you for your suggestion. I tried to implement the following:

transformed parameters {
    ...
    vector<lower=0, upper=1>[N] mu;
    vector<lower=1 / min(mu, 1 - mu)>[N] theta;
    ...
}

However, this results in a data compatibility error:

Ill-typed arguments supplied to function ‘min’:
(vector, vector)
Available signatures:
(int, int) => int
The first argument must be int but got vector
(array int) => int
Expected 1 arguments but found 2 arguments.
(array real) => real
Expected 1 arguments but found 2 arguments.
(vector) => real
Expected 1 arguments but found 2 arguments.
(row_vector) => real
Expected 1 arguments but found 2 arguments.

Any ideas to make overcome this? Would a custom “rowwise_min”-function work?

functions {
    vector bound(vector x) {
        vector[num_elements(x)] out;

        for (i in 1:num_elements(x)) {
            vector[2] v;
            v[1] = x[i];
            v[2] = 1 - x[i];

            out[i] = 1.0 / min(v);
        }

       return out;
    }
}

Or is there a better way to go?

7

Sorry for leading you astray! I should’ve realized that min wouldn’t work that way, because we don’t vectorize any binary functions that way!

You will need a function like you defined. The one you defined will work, but it’s inefficient to allocate and fill a vector, so I’d recommend just using

vector bound(vector x) {
  int N = cols(x);
  vector[N] y;
  for (n in 1:N) {
    y[n] = inv(min(x[n], 1 - x[n]));
  }
  return y;
}

You could use 1 / min(x[n], 1 - x[n]) rather than inv(min(x[n], 1 - x[n])), but inv is a bit more efficient.