2
$\begingroup$

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)
$\endgroup$
6
  • $\begingroup$ Perhaps this is useful? docs.pymc.io/Probability_Distributions.html $\endgroup$ Aug 26, 2021 at 16:11
  • $\begingroup$ Also, what is $m$ in your equation? I guess it's a typo for $n$. $\endgroup$ Aug 26, 2021 at 17:01
  • $\begingroup$ 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 $\endgroup$ Aug 26, 2021 at 17:05
  • 1
    $\begingroup$ Then $\frac{(1-\delta)}{2m}$ and $\frac{\delta}{2n}$ are misplaced and/or missing indices, surely? $\endgroup$ Aug 26, 2021 at 17:55
  • 2
    $\begingroup$ This is a programming question $\endgroup$
    – Firebug
    Aug 26, 2021 at 19:13

1 Answer 1

3
$\begingroup$

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?

$\endgroup$
3
  • $\begingroup$ 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) ? $\endgroup$ Aug 31, 2021 at 13:02
  • $\begingroup$ Good spot, I don't think that'll fix the derivate errors tho $\endgroup$ Aug 31, 2021 at 13:13
  • 1
    $\begingroup$ Hi @alanocallaghan Did you solve this? I need helpp on a similar question here: stats.stackexchange.com/questions/542138/… $\endgroup$
    – Tarik
    Sep 1, 2021 at 10:52

Not the answer you're looking for? Browse other questions tagged or ask your own question.