Search code examples
pymcpymc3

Model a distributions that is a combination of two others


I cannot find an example on this, or I don't see the similarity to my problem:

I tried to model a multimodal distribution that looks as if it were defined by the sum of two chi-squared distributions (though, one could take any combination of distributions ...).

Now I don't know how to "marry" those distributions. I am looking for something like this, but it doesn't work this way:

from pymc3 import Model, HalfNormal, find_MAP, sample, traceplot, ChiSquared, Deterministic basic_model = Model()

with basic_model:
    nu1 = HalfNormal("nu1", sd = 1)
    nu2 = HalfNormal("nu2", sd = 1)
    cs1 = ChiSquared("cs1", nu = nu1)
    cs2 = ChiSquared("cs2", nu = nu2)
    # this is wrong, but it shows what I would like to achieve:
    Y_obs = Deterministic("Y_obs", cs1, cs2, observed = tx)
    start = find_MAP(model = basic_model)
    trace = sample(2000, start = start)

traceplot(trace)

How would I do that? The resulting function should model something like this:

enter image description here


Solution

  • Check the following model, this model assumes your data is a mixture of two chi-squared distributions. The name for this type of model is a mixture model. Here, the categorical distribution is used to assign each data point to one of the two chi-squared distributions. You may need to adapt it to your problem, but I think is a good start.

    with pm.Model() as model:
        nus = pm.HalfCauchy('nus', beta=10, shape=2)
    
        category = pm.Categorical('category', p=[0.5, 0.5], shape=len(x))
    
        obs = pm.ChiSquared('obs', nu=nus[category], observed=x)
    
        step = pm.ElemwiseCategorical(vars=[category], values=[0, 1])
        trace = pm.sample(1000, step)