Combining two distributions

Questions
1

Hi there,

I though it was a simple matter but I’m getting nowhere with this. I’m basically “just” trying to construct an asymmetric StudentT distribution, i.e., two distributions with different sigmas, hoping to find a way to treat them as one and/or to normalize them so their peak value is the same. I hope that makes sense!

Anyway, using tt.switch works to choose the right distribution depending on whether the observed value is lower/greater than the model, like:
pm.Potential(label, tt.switch(tt.gt(observed, modeled), studentt_low, studentt_high))
but I’m not sure how the normalization would work here. Should I apply a factor to each?

An alternative would be to use a mixture, but I’m not sure how to combine HalfStudentT distributions or bounded StudentT distributions. Any help would be appreciated!

Cordially,
VL

1 Like
2

What is the “modeled” variable in your example doing?

3

It’s a tensor variable drawn from a grid of model data. But my question is general, how would you go about making an asymmetric Student-T distribution in PyMC3 with 2 different sigma and the same mean?

4

The reason I asked what the modeled is doing is because that might change how to approach your problem.

mu = pm.Normal('mu', 0, 1)
sigmas = pm.HalfNormal('sigmas', sigma=[1, 2], shape=2)  # You will probably need to use an ordered transform to avoid symmetries
sigma = tt.switch(tt.gt(observed, modeled), sigmas[0], sigmas[1])

obs = pm.StudentT('obs', mu=mu, sigma=sigma, observed=observed)

Were you looking for something similar in spirit to the example above? I am not sure if it will work out of the box, but maybe it gets you closer to your goal.

5

Yes, I was looking for something like that, I didn’t think of using the switch for the sigma, seems obvious now :slight_smile: I’ll make some tests and will let you know if there’s any issue. Thanks!

6

So a quick follow-up on that, as expected it returns an asymmetric Student-T, which is good, but it’s not normalized to the peak:

import pymc3 as pm
import theano.tensor as tt

modeled = 0

observed = -0.001
sigma = tt.switch(tt.gt(observed, modeled), 1., 20.)
print(pm.StudentT.dist(2, mu=modeled, sigma=sigma).logp(observed).eval())

observed = 0.001
sigma = tt.switch(tt.gt(observed, modeled), 1., 20.)
print(pm.StudentT.dist(2, mu=modeled, sigma=sigma).logp(observed).eval())

Returns:
-4.0354533
-1.0397215

I would assume the distribution has to be continuous, but maybe pymc3 does something under the hood I don’t understand.

7

It makes sense that the second one has a higher logp=-1.039... since it has lower standard deviation. The same way that a Normal distribution peak grows higher when you reduce the standard deviation.

If you really want the same logp at the peak when using different standard deviations I think you will have to derive a custom distribution that has a fixed height at zero and whose probability changes in the tails only (or do something funny with the degrees of freedom of the StudentT). This, however, would no longer be a StudentT, as far as I understand it.

8

You may also try the Asymmetric Laplace that was recently added to Pymc 3.11.0: Asymmetric Laplace distribution - Wikipedia. This distribution works by having two exponentials on either side, keeping the same peak but having different sigmas. Would that achieve what you need?

9

Hi Ricardo,

You’re right, at the end the distribution is not anymore a Student-T distribution, which was fine in my specific case but I should have been clearer.

I take the occasion to ask another question in this topic, I hope it’s ok. How would I go about drawing a Normal distribution but restricted to discrete values, for instance draw only integers for a Normal distribution with some mean and sigma?

Sorry if my questions are a bit random…!

Cheers,
Vian

10

May I ask with what goal in mind?

You can draw from a Normal distribution and round the values to integers, or you can use a discrete distribution with similar shape as the Normal (not so easy for small numbers). For example, both the Poisson with a very large mean or the Binomial with a large n and relatively mild p, approximate the Normal distribution.

No problem with random questions :slight_smile:

11

Thanks for your help. I’m basically trying to draw indices from a n-dimensional grid whose sides correspond to some physical parameters. Right now I draw each parameter as some kind of distribution, for instance Normal, and round the value, like you mentioned, either as integer or fractional value which I then use for nearest neighbor or linear interpolation. Sorry, I don’t have a simple code sample to extract.

It works fine but the trace then includes values which are actually not used to calculate the grid value at point. Of course I could always round the values in the trace as well a posteriori, but I was wondering whether it would just be simpler to draw directly these discrete values with a custom distribution.

edit: I could also use Deterministic to save the rounded values but that uses some useful memory.

12

You can save the rounded values in your model by wrapping them in a Deterministic:

indexes_continuous = pm.Normal(...)
indexes = pm.Deterministic('indexes', tt.round(indexes_continuous))

Then you can use indexes in your model and also check them in your trace without having to round a posteriori.