How to perform weighted inference

[mgregori]

I’m trying to do weighted inference for say a hierarchical model, using weight as proposed
here.
There’s a similar question here but the implementation the guy came up with it’s wrong in my opinion (instead of multiplying weights with likelihoods he is summing them).
Using a simplified model with regular inference like:

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

N,p=5000,2

X=np.random.normal(0,1,(N,p))
b=np.random.uniform(0,0.1,p)
y=np.random.poisson(np.exp(X.dot(b)),N)

w=np.random.uniform(0,1,N)

model = pm.Model()
with model:

    b0=pm.Normal('b0',mu=0,sd=.1)
    b1=pm.Normal('b1',mu=0,sd=.1)
    mu=pm.math.exp(b0*X[:,0]+b1*X[:,1])
    y_=pm.Poisson('Y_obs',mu=mu,observed=y)

    trace=pm.sample(1000)

I’ve tried to insert weights doing:

def weighted_logp(w,RV):
    def logp(x):
        return w*RV.logp(x)
    return logp

model = pm.Model()
with model:

    b0=pm.Normal('b0',mu=0,sd=.1)
    b1=pm.Normal('b1',mu=0,sd=.1)
    mu=pm.math.exp(b0*X[:,0]+b1*X[:,1])
    y_no_w=pm.Poisson('Y_obs',mu=mu,shape=N)
    y_w=pm.DensityDist('Y_w', weighted_logp(w,y_no_w), observed=y)

    trace=pm.sample(1000)

That is not working due to:

TypeError: can't turn [TensorConstant{[0. 0. 3. .. 2. 0. 3.]}] and {} into a dict. cannot convert dictionary update sequence element #0 to a sequence

With some tweaking, like:

def weighted_logp(w,RV):
    def logp(x):
        return w*RV.logp({RV.name:x})
    return logp

Still I couldn’t make it work.

ValueError: ('setting an array element with a sequence.', 'Container name "Y_obs"')

Any help?
I’ would like to make this as general as possible (i.e. avoid rewriting poisson likelihood)

Thank you!

[junpenglao]

If you want to scale the logp of the observed the same way as Bob’s solution in stackoverflow, the easiest way is to use a potential:

with model:

    b0=pm.Normal('b0',mu=0,sd=.1)
    b1=pm.Normal('b1',mu=0,sd=.1)
    mu=pm.math.exp(b0*X[:,0]+b1*X[:,1])
    y_=pm.Potential('Y_obs',w*pm.Poisson.dist(mu=mu).logp(y))

[mgregori]

Thank you for the prompt answer!

I’m not very familiar with the framework but Potential doesn’t allow to input observed values. Anyhow I’ve recycled you idea using DensityDist and I could make it work:

def weighted_logp(w,distributuion,**kwargs):
    return distributuion.dist(**kwargs).logp

model = pm.Model()
with model:

    b0=pm.Normal('b0',mu=0,sd=.1)
    b1=pm.Normal('b1',mu=0,sd=.1)
    mu=pm.math.exp(b0*X[:,0]+b1*X[:,1])
    y_w=pm.DensityDist('Y_w',weighted_logp(w,pm.Poisson,mu=mu) , observed=y)

    trace=pm.sample(1000)

Thanks again!

[junpenglao]

@mgregori

Potential doesn’t allow to input observed values

Just FYI, potential does not take observed values as input, but it can account for the observed value by adding the logp of the observation directly to the model log likelihood.

[gokl]

Hi,

@mgregori

There’s a similar question here but the implementation the guy came up with it’s wrong in my opinion (instead of multiplying weights with likelihoods he is summing them).

I think you are right. I added the normal Bernoulli likelihood with pm.Bernoulli and also added the weighted likelihood parts via pm.Potential. I should have used only the Potential likelihood. Thank you for pointing this out.

This would be the corrected example.

import numpy as np
import pymc3 as pm
import theano.tensor as T

N = 1000
P = 0.2

data = np.zeros(N)
data[:int(N * P)] = 1.0
weights  = np.ones(N)
weights[:int(N * P)] = 4.0

with pm.Model() as model:
    
    p_ = pm.Normal('p_',
                  mu=0.0,
                  sd=1000.0)
    p = pm.Deterministic('p', pm.math.invlogit(p_))
        
    pm.Potential('weights', 
             weights * pm.Bernoulli.dist(p=p).logp(data))

    trace = pm.sample(10000, cores=1, random_seed=1)
    
trace['p'].mean()  # 0.50019413

What do you think?

[mgregori]

Yup, I think that’s the right and easiest approach as junpenglao suggested.

On the other hand using DensityDist may allow you to then sample from the posterior distribution, unfortunately if you create wrapped functions and use multicores then some pickle issues may araise, like:

import numpy as np
import pymc3 as pm
import theano.tensor as tt
from functools import partial

N,p=5000,2

X=np.random.normal(0,1,(N,p))
b=np.random.uniform(0,0.1,p)
y=np.random.poisson(np.exp(X.dot(b)),N)

w=np.random.uniform(0,1,N)

def wlogp(x,w,distributuion,**kwargs):
return w*distributuion.dist(**kwargs).logp(x)
def weighted_distribution(w,distributuion,**kwargs):
return partial(wlogp,w=w,distributuion=distributuion,**kwargs),distributuion.dist(**kwargs).random

model = pm.Model()
with model:

b0=pm.Normal('b0',mu=0,sd=.1)
b1=pm.Normal('b1',mu=0,sd=.1)
mu=pm.math.exp(b0*X[:,0]+b1*X[:,1])
wlogp,wrandom=weighted_distribution(w,pm.Poisson,mu=mu)
y_w=pm.DensityDist('Y_w',wlogp, observed=y,random=wrandom)

trace=pm.sample(1000)

yields:

Auto-assigning NUTS sampler…
Initializing NUTS using jitter+adapt_diag…
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [b1, b0]
Could not pickle model, sampling singlethreaded.

But if I save the two helper functions in a separate package then multicore works. I wonder if makes sense to open an issue for that…

[byblian]

Is there any way to use the Potential to scale the likelihood with weights while also allowing imputation or handling of missing values in y?

[junpenglao]

Sure, just add one more potential term in the model.

[RusselDaries]

Hi @gokl , I tried sampling from posterior and I get back an empty dictionary. Is there a reason why? Your approach makes complete sense. Hence, my confusion if my implementation for posterior is incorrect.

Looking to you for advice.