Is this a censored-data problem or not?

Questions version agnostic
5

Thanks for your latest reply - I tried something else in the meantime:

Here’s my attempt at a custom likelihood for the censored data, in sample_censored_correct(). The rest is my toy problem for checking that naive modeling is biased.

"""Toy model for censored data"""

import numpy as np
import pymc3 as pm

true_rate = 0.25
N = 5_000

np.random.seed(45678)

customers = np.random.uniform(low=0, high=100, size=N)
complaints = np.random.poisson(lam=customers * true_rate)


def censor_data():
    measured = (customers >= 10) & ((complaints == 0) | (complaints >= 10))
    n = customers[measured]
    e = complaints[measured]
    return n, e


# noinspection PyTypeChecker, PyUnresolvedReferences
def do_sample_naive(n, e):
    assert len(n) == len(e)
    print(f"{len(n)} observations used")
    print(n[:20])
    print(e[:20])
    with pm.Model():
        r = pm.Exponential("r", 10)
        pm.Poisson("events", mu=n * r, observed=e)
        trace = pm.sample(2000, tune=1000)
        print(pm.summary(trace))


def sample_censored_wrong():
    n, e = censor_data()
    do_sample_naive(n, e)


def sample_uncensored():
    do_sample_naive(customers, complaints)


# noinspection PyTypeChecker,PyUnresolvedReferences
def sample_censored_correct():
    ns, es = censor_data()

    print(f"{len(ns)} observations used")
    print(ns[:20])
    print(es[:20])

    with pm.Model():
        r = pm.Exponential("r", 10)

        def truncated_poisson_like(e_, n_):
            lbda = n_ * r
            log_like = pm.Poisson.dist(mu=lbda).logp(e_)

            if e_ in [1, 2, 3]:
                return -np.inf  # log(0)
            else:
                # Normalization for the truncated values: subtract probabilities for 1, 2 and 3
                prob_sum = pm.math.exp(-lbda) * (lbda + lbda**2 / 2 + lbda**3 / 6)
                norm_constant = 1 - prob_sum
                return log_like - pm.math.log(norm_constant)

        pm.DensityDist(
            "observed", logp=truncated_poisson_like, observed={"e_": es, "n_": ns}
        )

        trace = pm.sample(2000, tune=1000)
        print(pm.summary(trace))


if __name__ == "__main__":
    print("\n\n🟢 Uncensored:\n")
    sample_uncensored()

    print("\n\n🔴 Censored, wrong:\n")
    sample_censored_wrong()

    print("\n\n🟣 Censored, correct:\n")
    sample_censored_correct()

The issue is that I get an error

Traceback (most recent call last):
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/pickle.py", line 603, in save
    self.save_reduce(obj=obj, *rv)
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/pickle.py", line 717, in save_reduce
    save(state)
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/pickle.py", line 560, in save
    f(self, obj)  # Call unbound method with explicit self
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/site-packages/dill/_dill.py", line 990, in save_module_dict
    StockPickler.save_dict(pickler, obj)
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/pickle.py", line 971, in save_dict
    self._batch_setitems(obj.items())
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/pickle.py", line 997, in _batch_setitems
    save(v)
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/pickle.py", line 603, in save
    self.save_reduce(obj=obj, *rv)
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/pickle.py", line 717, in save_reduce
    save(state)
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/pickle.py", line 560, in save
    f(self, obj)  # Call unbound method with explicit self
[... repeating ...]  
File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/pickle.py", line 560, in save
    f(self, obj)  # Call unbound method with explicit self
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/site-packages/dill/_dill.py", line 990, in save_module_dict
    StockPickler.save_dict(pickler, obj)
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/pickle.py", line 970, in save_dict
    self.memoize(obj)
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/pickle.py", line 510, in memoize
    self.write(self.put(idx))
  File "/Users/ehansis/anaconda/envs/pymc3/lib/python3.9/pickle.py", line 515, in put
    if self.proto >= 4:

RecursionError: maximum recursion depth exceeded in comparison

Any ideas on the error (or suggestions on my custom likelihood)? If not, I’ll look into your suggestion.

EDIT: I should mention that I’m on a quite old version of pymc3 (3.11.4), so maybe upgrading that would be a first step…