# Custom Likelihood Regression Pymc3 [closed]

I am struggling to implement a linear regression in pymc3 with a custom likelihood.

Suppose you have two independent variables $$x_1, x_2$$ and a target variable $$y$$, as well as an indicator variable $$\delta$$.

• When $$\delta$$ is 0, the likelihood function is standard least squares
• When $$\delta$$ is 1, the likelihood function is the least squares contribution only when the target variable is greater than the prediction

Mathematically, this is given as:

$$L(y_{pred}, y) = \frac{(1-\delta)}{2m}\sum_{i=1}^{m}(y_{pred_i} - y_i)^2 + \frac{\delta}{2n} \sum_{i=1}^{n} I(y_{pred_i}, y_i)(y_{pred_i} - y_i)^2$$

For indicator $$\delta$$, $$m$$ points with indicator $$\delta = 0$$, and $$n$$ points with indicator $$\delta = 1$$. Where $$I$$ is given by:

$$I(y_{pred_i}, y_i) = \begin{cases} 0 & y_i < y_{pred_i} \\ 1& y_i \geq y_{pred_i} \end{cases}$$

Example snippet of observed data:

x_1  x_2  š¯›æ   observed_target
10    1   0   100
20    2   0   50
5    -1   1   200
10   -2   1   100


Does anyone know how this can be implemented in pymc3? As a starting point...

model =  pm.Model()
with model as ttf_model:

intercept = pm.Normal('param_intercept', mu=0, sd=5)
beta_0 = pm.Normal('param_x1', mu=0, sd=5)
beta_1 = pm.Normal('param_x2', mu=0, sd=5)
std = pm.HalfNormal('param_std', beta = 0.5)

x_1 = pm.Data('var_x1', df['x1'])
x_2 = pm.Data('var_x2', df['x2'])

mu = (intercept + beta_0*x_0 + beta_1*x_1)

• Perhaps this is useful? docs.pymc.io/Probability_Distributions.html Aug 26, 2021 at 16:11
• Also, what is $m$ in your equation? I guess it's a typo for $n$. Aug 26, 2021 at 17:01
• Ah - m is the number of points you have with indicator (delta) value equal to 0, and n the number of points with delta equal to 1. Essentially average over each of the groups. So in the example data snipped I give - m would be 2 and n would be 2 also Aug 26, 2021 at 17:05
• Then $\frac{(1-\delta)}{2m}$ and $\frac{\delta}{2n}$ are misplaced and/or missing indices, surely? Aug 26, 2021 at 17:55
• This is a programming question Aug 26, 2021 at 19:13

This is as far as I've gotten:

import pymc3 as pm
import pandas as pd
import numpy as np
import theano.tensor as tt

with pm.Model() as model1:
intercept = pm.Normal('intercept', mu=0, sd=5)
beta_0 = pm.Normal('beta_0', mu=0, sd=5)
beta_1 = pm.Normal('beta_1', mu=0, sd=5)
std = pm.HalfNormal('std', sigma=0.5)
x_1 = pm.Data('x_1', data['x_1'])
x_2 = pm.Data('x_2', data['x_2'])
indicator = pm.Data('indicator', data['indicator'])
y_pred = (intercept + beta_0*x_1 + beta_1*x_2)
observed = pm.Normal("observed",
mu = y_pred, sigma = std,
observed = data['observed_target']
)
samples = pm.sample()


This runs, of course.

def logp(y_pred, y, indicator, total):
m = pm.math.sum(1 - indicator)
n = pm.math.sum(indicator)
greater = pm.math.gt(y, y_pred)
out = 0
out = pm.math.sum(
(
(indicator * greater) * pow(y_pred - y, 2)
/ (2 * pm.math.sum(indicator))
)
+ (
((1 - indicator) * pow(y_pred - y, 2))
/ (2 * pm.math.sum(1 - indicator))
)

)
return out

with pm.Model() as model2:
intercept = pm.Normal('intercept', mu=0, sd=5)
beta_0 = pm.Normal('beta_0', mu=0, sd = 5)
beta_1 = pm.Normal('beta_1', mu=0, sd = 5)
# std = pm.HalfNormal('std', sigma = 0.5)
x_1 = pm.Data('x_1', data['x_1'])
x_2 = pm.Data('x_2', data['x_2'])
indicator = pm.Data('indicator', np.array(data['indicator']))
y_pred = (intercept + beta_0*x_1 + beta_1*x_2)
observed = pm.DensityDist(
'observed',
logp,
observed = {
"y_pred": y_pred,
'y': data['observed_target'],
'indicator': indicator,
'total': len(data)
},
)
pm.sample()


This fails with errors:

ValueError: Mass matrix contains zeros on the diagonal.
The derivative of RV beta_1.ravel()[0] is zero.
The derivative of RV intercept.ravel()[0] is zero.


I might revisit. Otherwise, this might be a decent starting point?

• hey - thanks for this. Should this comparison be the other way around - greater = pm.math.gt(y_pred, y) ----> greater = pm.math.gt(y, y_pred) ? Aug 31, 2021 at 13:02
• Good spot, I don't think that'll fix the derivate errors tho Aug 31, 2021 at 13:13
• Hi @alanocallaghan Did you solve this? I need helpp on a similar question here: stats.stackexchange.com/questions/542138/… Sep 1, 2021 at 10:52