Censored pymc for power-law distribution of column density

Questions v5
, ,
11

The code for pymc censored is as follows:

def censored_regression(x, y, y_err):

    # data pre-processing
    left_censored = (y_err < -10)
    right_censored = (y_err > -10) & (y_err < -0)

    lower = [y[ind] if b else -np.inf for ind,b in enumerate(left_censored)]
    upper = [y[ind] if b else np.inf for ind,b in enumerate(right_censored)]

    with pm.Model() as model:

        A = pm.Normal('A', mu=0, sigma=10)
        B = pm.Normal('B', mu=0, sigma=10)
        σ = pm.HalfNormal("σ", sigma=1)

        normal_dist = pm.Normal.dist(mu=A * at.power(x, B), sigma=σ)
        pm.Censored("likelihood", normal_dist, lower=lower,
                    upper=upper, observed=y)

    with model:
        trace = pm.sample()

    return model, trace

The results are almost identical to the code above.

1 Like