How to implement bivariate poisson logbp

Questions v5
1

Hi, im trying to implement the bivariate poisson in pymc to use it in a CustomDist, but I’m struggling to work with TensorVariables because I have to compute the min over each vector input and then compute sum over a specific interval. I have tried to avoid this using .eval() but for lam1, lam2 and lam3 is not possible. Any suggestion? or hint to implement this direct or any idea?

The formula of the logp is:

\log(p(x, y) ) = -(\lambda_1+\lambda_2+\lambda_3) + x \log({\lambda_1}) - \log(x!) + y \log(\lambda_2) - \log(y!) + \\ \log \left( \sum_{i = 0}^{min(x,y)}{\binom xi \binom yi i! \left(\frac{\lambda_3}{\lambda_1\lambda_2}\right)^i} \right)

and the code I’ve been working on:

def bivpoiss_logp(value, lam1, lam2, lam3):
    x = value.eval()[0,:]
    y = value.eval()[1,:]

    logp = -(lam1 + lam2 + lam3) + x * np.log(lam1) - factln(x) + y * np.log(lam2) - factln(y)
    z = np.minimum(x,y)
    n = len(x)
    m = np.max(z)+1
    indexes = np.arange(m)
    term = np.array([i_term(x_i,y_i,indexes, lam1, lam2, lam3) for x_i,y_i in list(zip(x,y))])

    z = np.broadcast_to(z, (m,n))
    i_matrix = np.broadcast_to(np.arange(m), (n,m))
    valid_i = i_matrix <= z.T
    term[~valid_i] = 0
    result = np.sum(term, axis=1)
    return log(result).eval() + logp

and the i_term function define as follows:

def i_term(x,y, I, lambda_1, lambda_2, lambda_3):
    return comb(x, I)* comb(y, I) * factorial(I) * (lambda_3 / (lambda_1 * lambda_2)) ** I
2

You have to use PyTensor operations and you can’t use eval as that breaks the graph.

All np.foo should have an equivalent pt.foo, wherept. is an alias for pytensor.tensor

n = len(x) may need to be replaced by n = x.shape[0]

Your list comprehension can be avoided by making i_term vectorized.

term[~valid_i] = 0 has to be written as term = term[~valid_i].set(0)

I may have missed other hard bits.

1 Like
3

Thank you so much. I´ve alredy do it but I got another issue. The final logp looks like:

def bivpois_logp(value: TensorVariable, lam1: TensorVariable, lam2: TensorVariable, lam3: TensorVariable) -> TensorVariable:
    x = value[:,0]
    y = value[:,1]
    lam = lam3/(lam1*lam2)

    logp = -(lam1 + lam2 + lam3) + x * log(lam1) - factln(x) + y * log(lam2) - factln(y)
    z = pt.min(value, axis = 1)
    m = pt.max(z) +1 
    n = x.shape[0]
    indexes = pytensor.shared(np.arange(m.eval()))
    x_aux = x[:, None]
    y_aux = y[:, None]

    term = comb_pt(x_aux, indexes) * comb_pt(y_aux, indexes) * gamma(indexes+1) * ((lam)[:,None] ** indexes)

    z_aux = pt.zeros(shape = (n,m), dtype = "int") + z[:, None]
    ind_mat = pt.zeros(shape = (n,m), dtype = "int") + indexes
    valid_i = ind_mat <= z[:, None]
    term = term[~valid_i].set(0)
    result = pt.log(pt.sum(term, axis = 1))

    return logp + result

but only works when the variables are brodcastables, i mean, when value and lam shape is equal in one dimension, but running the next model shows that this is not always the case

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

    # Data
    home_team = pm.Data("home_team", home_idx, dims = "match")
    away_team = pm.Data("away_team", away_idx, dims = "match")

    # global model parameters
    home = pm.Normal("home", mu = 0, sigma = 1/0.0001)
    mu_att = pm.Normal('mu_att', mu = 0, sigma = 1/0.0001)
    mu_def = pm.Normal('mu_def', mu = 0, sigma = 1/0.0001)
    tau_att = pm.Gamma('tau_att', alpha = 0.1, beta = 0.1)
    tau_def = pm.Gamma('tau_def', alpha = 0.1, beta = 0.1)

    beta = pm.Normal('beta_con', mu = 0, sigma = 1/0.0001, dims = "match")

    # team model parameters
    att_star = pm.Normal("att_star", mu = mu_att, sigma = tau_att, dims = "team")
    def_star = pm.Normal("def_star", mu = mu_def, sigma = tau_def, dims = "team")

    # sum to zero
    att = pm.Deterministic("att", att_star - pm.math.mean(att_star), dims = "team")
    deff = pm.Deterministic("def", def_star - pm.math.mean(def_star), dims = "team")

    # Regression
    l1 = pt.exp(home + att[home_idx] + deff[away_idx])
    l2 = pt.exp(att[away_idx] + deff[home_idx])
    l3 = pt.exp(beta)

    # Likelihood
    BivariatePoisson = pm.CustomDist('goals_per_match', l1, l2, l3, logp = bivpois_logp, random = bivpois_rvs, dtype='int', observed = partidos_resultados[['home_goals', 'away_goals']], shape = (380, 2))

    trace = pm.sample(1000, tune=1500, cores = 4)

So my question would be: what shapes should I expect from l1,l2,l3 and the value for the logp?

4

This is odd, just do pt.arange(m)

The shapes of l1… l3 are the same as outside, you can check them by doing print(l1.eval().shape)

The shape of value is the same as the observed, which should be (380, 2)

1 Like
5

Ok thanks for the feedback. But I got this error

ValueError                                Traceback (most recent call last)
Cell In[106], line 2
      1 with model:
----> 2     trace = pm.sample(1000, tune=1500, cores = 4)
...
    456     )
    458 # inplace_pattern maps output idx -> input idx
    459 inplace_pattern = self.inplace_pattern

ValueError: Incompatible Elemwise input shapes [(1, 380), (1, 5)]

Could be the outputs’s shape of the function bivpoiss_logp but i can’t do .eval().shape inside of the logp function to check the shape because this error pops up

MissingInputError: Input 0 (def_star) of the graph (indices start from 0), used to compute Sum{axes=None}(def_star), was not provided and not given a value. Use the PyTensor flag exception_verbosity='high', for more information on this error.

Can you give me any advice? please

Additionally I run model.debug(verbose = True) and I got this:

point={'home': array(0.), 'mu_att': array(0.), 'mu_def': array(0.), 'tau_att_log__': array(0.), 'tau_def_log__': array(0.), 'beta_con': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0.]), 'att_star': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0.]), 'def_star': array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0.])}

The variable goals_per_match has the following parameters:
0: Exp [id A] <Matrix(float64, shape=(1, ?))>
 └─ Add [id B] <Matrix(float64, shape=(1, ?))>
    ├─ ExpandDims{axes=[0, 1]} [id C] <Matrix(float64, shape=(1, 1))>
    │  └─ home [id D] <Scalar(float64, shape=())>
    ├─ ExpandDims{axis=0} [id E] <Matrix(float64, shape=(1, ?))>
    │  └─ AdvancedSubtensor1 [id F] <Vector(float64, shape=(380,))>
    │     ├─ Sub [id G] <Vector(float64, shape=(?,))> 'att'
    │     │  ├─ att_star [id H] <Vector(float64, shape=(?,))>
    │     │  └─ True_div [id I] <Vector(float64, shape=(1,))>
    │     │     ├─ ExpandDims{axis=0} [id J] <Vector(float64, shape=(1,))>
    │     │     │  └─ Sum{axes=None} [id K] <Scalar(float64, shape=())>
    │     │     │     └─ att_star [id H] <Vector(float64, shape=(?,))>
    │     │     └─ Cast{float64} [id L] <Vector(float64, shape=(1,))>
    │     │        └─ MakeVector{dtype='int64'} [id M] <Vector(int64, shape=(1,))>
    │     │           └─ Shape_i{0} [id N] <Scalar(int64, shape=())>
    │     │              └─ att_star [id H] <Vector(float64, shape=(?,))>
    │     └─ [ 7  2  8 ... 5 13 11] [id O] <Vector(uint8, shape=(380,))>
    └─ ExpandDims{axis=0} [id P] <Matrix(float64, shape=(1, ?))>
       └─ AdvancedSubtensor1 [id Q] <Vector(float64, shape=(380,))>
          ├─ Sub [id R] <Vector(float64, shape=(?,))> 'def'
          │  ├─ def_star [id S] <Vector(float64, shape=(?,))>
          │  └─ True_div [id T] <Vector(float64, shape=(1,))>
          │     ├─ ExpandDims{axis=0} [id U] <Vector(float64, shape=(1,))>
          │     │  └─ Sum{axes=None} [id V] <Scalar(float64, shape=())>
          │     │     └─ def_star [id S] <Vector(float64, shape=(?,))>
          │     └─ Cast{float64} [id W] <Vector(float64, shape=(1,))>
          │        └─ MakeVector{dtype='int64'} [id X] <Vector(int64, shape=(1,))>
          │           └─ Shape_i{0} [id Y] <Scalar(int64, shape=())>
          │              └─ def_star [id S] <Vector(float64, shape=(?,))>
          └─ [ 0  1  6 ... 18  9 10] [id Z] <Vector(uint8, shape=(380,))>
1: Exp [id BA] <Matrix(float64, shape=(1, ?))>
 └─ Add [id BB] <Matrix(float64, shape=(1, ?))>
    ├─ ExpandDims{axis=0} [id BC] <Matrix(float64, shape=(1, ?))>
    │  └─ AdvancedSubtensor1 [id BD] <Vector(float64, shape=(380,))>
    │     ├─ Sub [id G] <Vector(float64, shape=(?,))> 'att'
    │     │  └─ ···
    │     └─ [ 0  1  6 ... 18  9 10] [id Z] <Vector(uint8, shape=(380,))>
    └─ ExpandDims{axis=0} [id BE] <Matrix(float64, shape=(1, ?))>
       └─ AdvancedSubtensor1 [id BF] <Vector(float64, shape=(380,))>
          ├─ Sub [id R] <Vector(float64, shape=(?,))> 'def'
          │  └─ ···
          └─ [ 7  2  8 ... 5 13 11] [id O] <Vector(uint8, shape=(380,))>
2: Exp [id BG] <Matrix(float64, shape=(1, ?))>
 └─ ExpandDims{axis=0} [id BH] <Matrix(float64, shape=(1, ?))>
    └─ beta_con [id BI] <Vector(float64, shape=(?,))>
The parameters evaluate to:
0: [[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]
1: [[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]
2: [[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]
The variable goals_per_match logp method raised the following exception: Input dimension mismatch: (input[%i].shape[%i] = %lld, input[%i].shape[%i] = %lld)

Input dimension mismatch: (input[%i].shape[%i] = %lld, input[%i].shape[%i] = %lld)
Apply node that caused the error: Composite{((i2 / (i0 * i1)) ** i3)}(ExpandDims{axis=1}.0, ExpandDims{axis=1}.0, ExpandDims{axis=1}.0, [[[0 1 2 3 4]]])
Toposort index: 22
Inputs types: [TensorType(float64, shape=(1, 1, None)), TensorType(float64, shape=(1, 1, None)), TensorType(float64, shape=(1, 1, None)), TensorType(int64, shape=(1, 1, 5))]
Inputs shapes: [(1, 1, 380), (1, 1, 380), (1, 1, 380), (1, 1, 5)]
Inputs strides: [(3040, 3040, 8), (3040, 3040, 8), (3040, 3040, 8), (40, 40, 8)]
Inputs values: ['not shown', 'not shown', 'not shown', array([[[0, 1, 2, 3, 4]]])]
Outputs clients: [[Composite{((i0 * i1 * i2) / i3)}([[[1.0000e ... 000e+01]]], [[[  2.] ... [ 24.]]], Composite{((i2 / (i0 * i1)) ** i3)}.0, [[[  2.  i ... 6.  24.]]])]]

HINT: Re-running with most PyTensor optimizations disabled could provide a back-trace showing when this node was created. This can be done by setting the PyTensor flag 'optimizer=fast_compile'. If that does not work, PyTensor optimizations can be disabled with 'optimizer=None'.
HINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.

I think I know what’s going on, when I implented the vectorize form I was thinking of numpy arrays so when I compute lam ** indexes with lam’s shape (1,n) and indexes’s shape (m,1) I’ve end up with lam ** indexes 's shape (m,n) and like the i,j element is lam[i] ** indexes[j]. It is possible with TensorVariables?

6

Yes, works the same with TensorVariables

7

Ok thanks, I already checked like this

index = pt.arange(10)

lams = pm.Exponential.dist(lam = 0.5, shape = (5,3))
lam1 = lams[:,0]
lam2 = lams[:,1]
lam3 = lams[:,2]
lam = (lam3 / (lam2 * lam1))[:,None]

resultado = lam ** index

And the print looks nice:

[0 1 2 3 4 5 6 7 8 9] (10,)
[[  3.92221234]
 [  6.99686109]
 [  0.72190472]
 [  1.06022012]
 [695.87921068]] (5,1)
[[1.00000000e+00 3.92221234e+00 1.53837496e+01 6.03383326e+01
  2.36659753e+02 9.28229802e+02 3.64071438e+03 1.42796549e+04
  5.60078385e+04 2.19674635e+05]
 [1.00000000e+00 6.99686109e+00 4.89560650e+01 3.42538786e+02
  2.39669631e+03 1.67693511e+04 1.17332820e+05 8.20961444e+05
  5.74415318e+06 4.01910418e+07]
 [1.00000000e+00 7.21904723e-01 5.21146429e-01 3.76218069e-01
  2.71593601e-01 1.96064703e-01 1.41540035e-01 1.02178420e-01
  7.37630839e-02 5.32499186e-02]
 [1.00000000e+00 1.06022012e+00 1.12406671e+00 1.19175814e+00
  1.26352596e+00 1.33961565e+00 1.42028747e+00 1.50581735e+00
  1.59649785e+00 1.69263915e+00]
 [1.00000000e+00 6.95879211e+02 4.84247876e+05 3.36978030e+08
  2.34496005e+11 1.63180895e+14 1.13554192e+17 7.90200018e+19
  5.49883765e+22 3.82652680e+25]] (5, 10)

But why I keep getting this error when I run the model and the logbp looking the same as the example above? And also why I can’t check the value of lam1, lam2 and lam3 inside of the logp using the function .eval()?

8

Because the logp graph is made of the values of the RVs, so you would need to propose a value for them when you want to evaluate an expression inside the logp function. However this is not actually allowed right now, I opened this pull-request to allow that in the future: Allow debug evaling IR logp graphs by ricardoV94 · Pull Request #7666 · pymc-devs/pymc · GitHub

1 Like
9

In the meantime, if you can provide a fully reproducible snippet with your latest attempts I may be able to help.

1 Like
10

Wow thank you so much, it will be super useful.

Yes of course, I made it a google colab: Google Colab

11

Your lam shapes don’t make sense?

image

That would fail as soon as you do lam3/(lam1*lam2) in your logp (just like they would in numpy)

12

Besides the shape mismatch, you probably will want to convert many of your expressions to log scale, for stability.

All this using, using addition / subtraction with gammaln, instead of gamma (computing the output on log scale)

def comb_pt(k: TensorVariable, n:TensorVariable) -> TensorVariable:
    res = gamma(k+1) / (gamma(k-n+1) * gamma(n +1))
    return res

Power to multiplication of logs, multiplication to addition of logs…

    pre2 = (lam)[:, None] ** indexes
    pre1 = comb_pt(x_aux, indexes) * comb_pt(y_aux, indexes) * gamma(indexes+1)
    ...
    term = pre2*pre1
13

Sorry I forgot to undo the comment when I use index the prior to get the shape of (380,), but now it’s all good and the error still if you want to check it now.

and then like computnig the exp of it? to keep the correct values of the logp I mean

14

Well the logp is on log-scale, so probably you can do without any exps in it. May need to use stuff like logaddexp, logdiffexp, logsumexp

1 Like
15

I’ll check the NB later. If we manage to get it working would be great to add it to pymc-extras, if you are so inclined

1 Like
16

Realized a simple shortcut to testing the logp and shapes, is to call the logp function yourself with the model variables :wink:

bivpois_logp(BivariatePoisson, l1, l2, l3)

It will use the random variables, not their values, but should be good to check shape errors.

1 Like
17

I fixed a couple of things in the notebook.

Namely providing a signature to CustomDist to let it know it’s a multivariate distribution, and using dist instead of random, since you were using PyMC variables, not numpy methods:

    BivariatePoisson = pm.CustomDist(
        'goals_per_match', 
        l1, l2, l3, 
        logp = bivpois_logp, 
        dist = bivpois_rvs, 
        signature="(),(),()->(2)", 
        observed = games[['score_home', 'score_away']],
    )

Also you’re not allowed to return multiple variables in the dist/random function so I changed it into:

def bivpois_rvs(lam1, lam2, lam3, size=None):
    z1 = pm.Poisson.dist(mu = lam1, size = size)
    z2 = pm.Poisson.dist(mu = lam2, size = size)
    z3 = pm.Poisson.dist(mu = lam3, size = size)

    x1 = z1 + z3
    x2 = z2 + z3

    return pt.stack([x1, x2], axis=-1)

Now there are no shape errors, but the model has a nan logp :slight_smile:

image

Updated notebook: Google Colab

1 Like
18

Got a version of the logp to not output nan, but did not check carefully wether it is correct:

from pymc.math import log
from pymc.distributions.dist_math import factln, logpow
import pytensor.tensor as pt
from pytensor.tensor import TensorVariable
from pytensor.tensor.math import gammaln

def comb_pt(k: TensorVariable, n:TensorVariable) -> TensorVariable:
    # res = gamma(k+1) / (gamma(k-n+1) * gamma(n +1))
    res = gammaln(k+1)- (gammaln(k-n+1) + gammaln(n+1))
    return res
    
def bivpois_logp(value: TensorVariable, lam1: TensorVariable, lam2: TensorVariable, lam3: TensorVariable) -> TensorVariable:
    x = value[:,0]
    y = value[:,1]
    lam = lam3/(lam1*lam2)

    first_term = -(lam1 + lam2 + lam3) + x * log(lam1) - factln(x) + y * log(lam2) - factln(y)

    # Sum terms
    z = pt.min(value, axis = -1)

    # Do square sum on the largest z
    m = pt.max(z) + 1
    indexes = pt.arange(m)

    sum_terms = (
        comb_pt(x[:, None], indexes) 
        + comb_pt(y[:, None], indexes) 
        + gammaln(indexes+1)
        + logpow(lam[:,None], indexes)
    )

    # Set to zero sum terms that go beyond z
    sum_terms = sum_terms[z[:, None] > indexes].set(-np.inf)

    logp = first_term + pt.logsumexp(sum_terms, axis = -1)
    return logp

It’s in the same Colab copy I linked above.

Sampling seems to go at snail pace though. Not sure if the model is well specified (and much less sure I implemented the logp correctly, do triple-check it!)

1 Like
19

Thank you so much for your help, the final version of the logp is the next one:

def comb_pt(k: TensorVariable, n:TensorVariable) -> TensorVariable:
    return gammaln(k+1) - gammaln(k-n+1) - gammaln(n+1)

def bivpois_logp(value: TensorVariable, lam1: TensorVariable, lam2: TensorVariable, lam3: TensorVariable) -> TensorVariable:
    x = value[:,0]
    y = value[:,1]
    lam = lam3/(lam1*lam2)

    first_term = - (lam1 + lam2 + lam3) + x * log(lam1) - factln(x) + y * log(lam2) - factln(y)

    z = pt.min(value, axis = -1)
    m = pt.max(z) + 1
    indexes = pt.arange(m)

    sum_terms = (
        comb_pt(x[:, None], indexes)
        + comb_pt(y[:, None], indexes)
        + gammaln(indexes+1)
        + log(lam)[:,None] * indexes
    )

    sum_terms = sum_terms[z[:, None] > indexes].set(-np.inf)

    logp = first_term + pt.logsumexp(sum_terms, axis = -1)
    return logp

Yes, would be awesome, How can I help with it?

20

Are you sure about this? Aren’t the first z terms what you want to keep? You are erasing them this way