Combining two distributions

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(, 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!


1 Like

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

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?

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(, 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.

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!

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(, modeled), 1., 20.)
print(pm.StudentT.dist(2, mu=modeled, sigma=sigma).logp(observed).eval())

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


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

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.

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?

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…!


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:

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.

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.