Censored distributions

[domenzain]

Is there a way to express left, right or interval censored distributions?

Here's a small example illustrating the differences with numpy:

#!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(123)

# Create normally distributed samples.
size = 10000
sigma = 1.
mu = 1.
samples = np.random.normal(mu, sigma, size)

# Define bounds for the distribution
lower=0.
upper=1.

# Keep information about samples being outside of the range.
censored = samples.copy()
censored[censored > upper] = upper
censored[censored < lower] = lower

# Lose all information outside the range that can be sampled. 
truncated = samples[(samples >= lower) & (samples <= upper)]

plt.hist([samples, censored, truncated], bins=40)
plt.show()

Using Bound and Normal I can describe something like the truncated distribution.
But many processes behave like censored. How can the latter be expressed in PyMC3?

[fonnesbeck]

We do not have built-in support for censoring. When I have to deal with this (for survival analysis), I typically just use a vector of indicators to indicate censoring.

[aseyboldt]

There is a section about this in the stan manual (section 11.3, Cencored data), and this should work pretty much the same in pymc3. For the second method you might have to code the ccdf functions on your own (a lot of special functions are in theano), as far as I know pymc does not provide those out of the box. You can use a Potential to do the same as target += something.

[domenzain]

I'll give it a go and create an example.

[domenzain]

The first technique in the manual was straightforward to test. And seems to work as expected, despite being very heavy computationally. I'll give the second technique a go, too.

I've also found a paper Likelihood Estimation for Censored
Random Vectors (PDF, Schnedler, 2005) that describes a technique to automatically generate the likelihood of a broad category of censoring problems. Perhaps this could be implemented eventually?

When comparing with a Normal model and uncensored data, and Bounded Normal with truncated data, the latter behaves differently from what I'd expect. The estimate of both mu and sigma is off by a factor of 2.

Am I doing something obviously incorrect here? Did I forget to normalize something?

#!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm

# Produce normally distributed samples
np.random.seed(123)
size = 10000
sigma = 5.
mu = 13.
samples = np.random.normal(mu, sigma, size)
expected_mu = np.mean(samples)
expected_sigma = np.std(samples)

# Set censoring/truncation limits.
high = 10.
low = 0.

# Truncate samples
truncated = samples[(samples > low) & (samples < high)]

# Keep tabs on left/right censoring for censored model
n_right_censored = len(samples[samples >= high])
n_left_censored = len(samples[samples <= low])
n_observed = len(samples) - n_right_censored - n_left_censored

# Omniscient model
with pm.Model():
    mu = pm.Normal('mu', mu=0., sd=(high-low)/2.)
    sigma = pm.HalfNormal('sigma', sd=(high-low)/2.)
    observed = pm.Normal('observed', mu=mu, sd=sigma, observed=samples)
    trace = pm.sample(10000, njobs=6)
    print('Omniscient model')
    pm.summary(trace)
    pm.traceplot(trace)
    plt.show()

# Truncated model
with pm.Model():
    mu = pm.Normal('mu', mu=0., sd=(high-low)/2.)
    sigma = pm.HalfNormal('sigma', sd=(high-low)/2.)
    observed = pm.Bound(
        pm.Normal, lower=low, upper=high
    )('observed', mu=mu, sd=sigma, observed=truncated)
    trace = pm.sample(10000, njobs=6)
    print('Truncated model')
    pm.summary(trace)
    pm.traceplot(trace)
    plt.show()

# Censored model
with pm.Model():
    mu = pm.Normal('mu', mu=0., sd=(high-low)/2.)
    sigma = pm.HalfNormal('sigma', sd=(high-low)/2.)
    right_censored = pm.Bound(pm.Normal, lower=high)(
        'right_censored', mu=mu, sd=sigma, shape=n_right_censored
    )
    left_censored = pm.Bound(pm.Normal, upper=low)(
        'left_censored', mu=mu, sd=sigma, shape=n_left_censored
    )
    observed = pm.Normal(
        'observed',
        mu=mu,
        sd=sigma,
        observed=truncated,
        shape=n_observed
    )
    trace = pm.sample(300)
    print('Censored model')
    pm.summary(trace)
    pm.traceplot(trace)
    plt.show()

[aseyboldt]

Bound does not do what you think it does. It does not correct the logp with the cdf or ccdf of the distribution, but only takes care of re re-parametrization. I'd consider this a bug in pymc3, but others might disagree. At the moment you can work around that by manually adding a potential for the cdf (I think something likepm.Potential('normal_truncation', - tt.log(tt.erf((high - mu) / sigma) - tt.erf((low - mu) / sigma))), but I haven't tested it)

[aseyboldt]

This seems to work:

# Truncated model
with pm.Model() as model:
    mu = pm.Flat('mu')
    sigma = pm.HalfCauchy('sigma', beta=2.5)
    pm.Normal('observed', mu=mu, sd=sigma, observed=truncated)
    trunc = tt.log(tt.erf((high - mu) / sigma / np.sqrt(2)) - tt.erf((low - mu) / sigma / np.sqrt(2)))
    pm.Potential('normal_truncation', - len(truncated) * trunc)

    trace = pm.sample(2000, tune=1000, njobs=4, init=None)

[domenzain]

Thank you for that.

Here's the second technique in the manual for censored data. Unsurprisingly this is much better computationally.

It would have been much quicker had the standard functions been around already (the numerical stability is not great for lccdf using the naive implementation). If anything providing the functions needed to fix Bound in #1843 would be a huge convenience.

import theano.tensor as tt

# Unimputed censored model
def normal_lcdf(mu, sigma, x):
    z = (x - mu) / tt.sqrt(2. * sigma ** 2.)
    return tt.switch(
        tt.lt(z, -1.0),
        tt.log(tt.erfcx(-z/tt.sqrt(2.))/2.) - tt.sqr(tt.abs_(z))/2,
        tt.log1p(-tt.erfc(z/tt.sqrt(2.))/2.)
    )


def normal_lccdf(mu, sigma, x):
    z = (x - mu) / tt.sqrt(2. * sigma ** 2.)
    return tt.switch(
        tt.gt(z, 1.0),
        tt.log(tt.erfcx(z/tt.sqrt(2.))/2) - tt.sqr(tt.abs_(z))/2.,
        tt.log1p(-tt.erfc(-z / tt.sqrt(2.))/2.)
    )


def censored_right_likelihood(mu, sigma, n_right_censored, upper_bound):
    return n_right_censored * normal_lccdf(mu, sigma, upper_bound)


def censored_left_likelihood(mu, sigma, n_left_censored, lower_bound):
    return n_left_censored * normal_lcdf(mu, sigma, lower_bound)

with pm.Model():
    mu = pm.Normal('mu', mu=0., sd=(high-low)/2.)
    sigma = pm.HalfNormal('sigma', sd=(high-low)/2.)
    observed = pm.Normal(
        'observed',
        mu=mu,
        sd=sigma,
        observed=truncated,
    )
    right_censored = pm.Potential(
        'right_censored',
        censored_right_likelihood(mu, sigma, n_right_censored, high)
        )
    left_censored = pm.Potential(
        'left_censored',
        censored_left_likelihood(mu, sigma, n_left_censored, low)
    )
    trace = pm.sample(50000)
    print('Censored model')
    pm.summary(trace)
    pm.traceplot(trace)
    plt.show()

[domenzain]

Should I make a pull request for an example on the different methods of dealing with censored data?

[aseyboldt]

I guess the pull request is a question for @fonnesbeck or @twiecki, but I doubt they'd say no. But it might be better to wait for #1843.
Good idea about the numerical stability of normal_cdf, however the if will not do what you want. tt.lt and tt.gt return variables, that represent the comparison in a computational graph, not booleans. So python will always take the first branch, as the return value of tt.gt does not evaluate to False or None. You can use tt.switch to do what you want.

[domenzain]

Thank you again. I've updated the code above to use tt.switch.

[twiecki]

Yes, I think such a PR would be great. It's a bit unfortunate we have to add CDFs too though.

[domenzain]

I'll make a PR with these examples for PyMC3 as it stands, with the log CDF and log CCDF hard coded in the example as above.

Where would you put these functions otherwise? As methods of each distribution?

[twiecki]

Right, exactly.

…

On Fri, Mar 3, 2017 at 2:47 PM, domenzain ***@***.***> wrote: I'll make a PR with these examples for PyMC3 as it stands, with the log CDF and log CCDF hard coded in the example as above. Where would you put these functions otherwise? As methods of each distribution? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1833 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AApJmBabUoSvr4uF65itI-DqTyd-GLogks5riBnxgaJpZM4MNTZm> .

[fonnesbeck]

It might be handy to have log-CDFs for our distributions. Definitely not a bad thing.

As far as an interface goes, it might be nice to have it mirror Bounded, such as:

censored_gamma = pm.Censored(pm.Gamma, right=censoring_times)