Custom Likelihood Regression Pymc3 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)


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

data = pd.read_csv("pymc.csv")


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?
