# How to perform weighted inference

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!

2 Likes

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))
``````
1 Like

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!

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.

2 Likes

Hi,

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?

4 Likes

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…