When I use the following code (my priors and all are a bit different so the numbers are bigger)

```
def gamma_noise_regression(mu, alpha, beta, size=None):
return mu + pm.Gamma.dist(alpha=alpha, beta=beta, size=size)
```

The I get the following distribution of values when making predictions on new data:

The vertical line passes at 4020, which is where my truncation point is (values below this threshold are infeasible)

However, when I use the following code (capping `mu`

at 4020):

```
def gamma_noise_regression(mu, alpha, beta, size=None):
return pm.math.maximum(mu, 4020) + pm.Gamma.dist(alpha=alpha, beta=beta, size=size)
```

I get the following distribution of predicted values (along with about 2% divergences):

As you can see, the minimum of the values is quite far from the truncation point. What I want is to get the same distribution as in the first figure, but cut off at the truncation point (4020).

Also, when I use the following code:

```
def gamma_noise_regression(mu, alpha, beta, size=None):
return 4020 + pm.math.log1pexp(mu) + pm.Gamma.dist(alpha=alpha, beta=beta, size=size)
```

to avoid a hard boundary, I get sampling errors.

Here is my full code

```
def gamma_noise_regression(mu, alpha, beta, size=None):
return 4020 + pm.math.log1pexp(mu) + pm.Gamma.dist(alpha=alpha, beta=beta, size=size)
with pm.Model() as model:
# Define priors
const = pm.Normal("const",
mu=7_000,
sigma=1_000)
coef = pm.Uniform("coef", lower=-1, upper=0.0)
mean = X * coef + const
gamma_alpha = pm.Gamma("alpha", alpha=200.0, beta=10.0)
gamma_beta = pm.Gamma("beta", alpha=100.0, beta=10.0)
# Define likelihood
likelihood = pm.CustomDist('y', mean, gamma_alpha, gamma_beta,
dist=gamma_noise_regression,
observed=y)
# Inference!
# draw posterior samples using NUTS sampling
chain_number = 2
idata = pm.sample(SAMPLE_NUMBER//chain_number,
chains=chain_number,
tune=5_000,
)
# SamplingError: Initial evaluation of model at starting point failed!
# Starting values:
# {'const': array(7000.75080922), 'coef_interval__': array(-0.98534689), 'alpha_log__': array(3.41875223), 'beta_log__': array(2.06907241)}
# Logp initial evaluation results:
# {'const': -7.83, 'coef': -1.62, 'alpha': -18.98, 'beta': -1.14, 'y': -inf}
```

EDIT: I don’t have the proper vocabulary to explain this, but the issue with the results in the second figure is that in a case where `mu`

is sampled to be 1,000 and the Gamma distribution returns a sample of 4,000, the sum `mu+gamma`

is feasible under my model, but capping `mu`

at a value of 4,020 makes that sum infeasible (since `mu`

is smaller than 4,020 in this example). It is the sum that needs to be capped. That is why I feel implementing my own `logp`

for this sum is my best bet, but I don’t know how to do that