I am now trying to fit a regression model with my own function such as:
def fn(a,b,c,x1, x2):
y = 1/a + a/bx1 + c/ax2
y should always be positive, and all my predictors, x1 and x2 are positive. The parameters, a,b, and c can be negative, and
my observed values (likelihood functions) follow the gamma or wald distributions, but as you know that these distributions cannot have negative values.
To fit this model, I tried to do something like this:
with pm.Model() as model:
a = pm.Normal(‘a’, mu=0, sigma=100)
b = pm.Normal(‘a’, mu=0, sigma=100)
c = pm.Normal(‘a’, mu=0, sigma=100)
mu = fn(a,b,c,x1,x2)
ϵ = pm.HalfCauchy(‘ϵ’, 1)
y_pred = pm.Gamma(‘y_pred’, mu=mu, sigma=ϵ, observed=y_obs)
trace = pm.sample(100000, tune=50000, chains=2, target_accept = 0.999)
With this setting, I got an error because my function, fn, can produce a negative value, but mu cannot be negative in the Gamma likelihood function. In this case, how can I regulaize or penalize the model? or what is a better way to achieve this task?
Thank you for your time and help.