Multiplicative Factor Model (degradations)

Questions version agnostic
2

I found the basic algorithm elsewhere on discourse, and had to update it a bit. So now I have a shifted Exponential distribution, as described here: Exponential likelihood at different location. Hurray.

But I really want to flip the distribution, so I define this version, and I can sample from it.

import pymc as pm
from pytensor.tensor import TensorVariable

def dist(lam: TensorVariable, shift: TensorVariable, size: TensorVariable) -> TensorVariable:
    return -pm.Exponential.dist(lam, size=size) + shift

with pm.Model() as m:
    lam = pm.HalfNormal("lam")
    shift = -1
    pm.CustomDist(
        "custom_dist",
        lam, shift, dist=dist,
        # observed=[-1, -1, 0],
        observed=[-2, -2, -1], # Make sure observed data is valid (not 0 probability)
    )

    prior = pm.sample_prior_predictive()
    posterior = pm.sample()

for shift in [-1, 2, 3]:
  x = pm.CustomDist.dist(4, shift, dist=dist)
  x_draws = pm.draw(x, draws=4000)
  plt.hist(x_draws, bins=20, label=shift)
plt.legend();

shifted_exp

All is good.

My real data looks like this, adding the linear regression first.

size = 200
true_intercept = 4
true_slope = 6

x = np.linspace(0, 1, size)
# y = a + b*x
true_regression_line = true_intercept + true_slope * x
# add noise, rate = 1/scale
y = true_regression_line * 1/(1+np.random.exponential(scale=.5, size=size))

# Here's what it looks like
plt.figure(figsize=(3,3))
plt.scatter(x, y)
plt.plot(x, true_regression_line, color='red', linewidth=3)
plt.xlabel('x'); plt.ylabel('y');

sample_data

But when I try to do inference it says: RuntimeError: The logprob terms of the following value variables could not be derived:

Here is the code I’m running. Why can’t it derive the logp? And how do I do it myself?

with pm.Model() as linreg3:
  lam = pm.HalfNormal('sigma', initval=0.0) # , beta=10, initval=1.)
  intercept = pm.Normal('Intercept', 0, sigma=20)
  x_coeff = pm.Normal('Slope', 0, sigma=20)

  likelihood = pm.CustomDist('y', lam, 
                              intercept + x_coeff * x, 
                              dist=dist,
                              observed=y, 
                              )
  trace2 = pm.sample()

I’m confused about what I broke, and I’m not sure how to proceed.

Thanks for any help.

– Malcolm