Data aggregation in crowdsensing application

Questions
18

@junpenglao I’m almost finishing my little project, in the last version I’m trying to use in a loop the last posterior as the new prior, in particular I want to update the prior for the standard deviation. This is the code:

import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
from pymc3.distributions import Interpolated
from scipy import stats

def from_posterior(param, samples):
    smin, smax = np.min(samples), np.max(samples)
    width = smax - smin
    x = np.linspace(smin, smax, 100)
    y = stats.gaussian_kde(samples)(x)

    x = np.concatenate([[x[0] - 3 * width], x, [x[-1] + 3 * width]])
    y = np.concatenate([[0], y, [0]])
    return Interpolated(param, x, y)

# Initialize random number generator
np.random.seed(123)

# Target values
t = [38.095, 13.346, 30.106, 13.353]

# Users' precision in term of sigma
s = [0.001, 0.005, 0.003, 0.012]

# Number of measures for each target
n_measures = 10

# Simulate some data
measures = []
for i in range(0,4):
    y1 = t[i] + np.random.randn(n_measures) * s[0]
    y2 = t[i] + np.random.randn(n_measures) * s[1]
    y3 = t[i] + np.random.randn(n_measures) * s[2] 
    y4 = t[i] + np.random.randn(n_measures) * s[3]
    x = np.array([y1, y2, y3, y4])
    measures.append(x) 

with pm.Model() as hierachical_model:
    # Latent mean and sigma, I want to estimate these parameters
    mu = pm.Normal('mu', 0., 10.)
    rho = pm.HalfNormal('rho', 10.)

    sigma1 = pm.HalfNormal('sigma1', rho)
    sigma2 = pm.HalfNormal('sigma2', rho)
    sigma3 = pm.HalfNormal('sigma3', rho)
    sigma4 = pm.HalfNormal('sigma4', rho)

    like1 = pm.Normal('like1', mu, sigma1, observed=measures[0][0])
    like2 = pm.Normal('like2', mu, sigma2, observed=measures[0][1])
    like3 = pm.Normal('like3', mu, sigma3, observed=measures[0][2])
    like4 = pm.Normal('like4', mu, sigma4, observed=measures[0][3])

    trace = pm.sample(1000)

for j in range(1,4):

    with pm.Model() as model:
        # Latent mean and sigma, I want to estimate these parameters
        mu = pm.Normal('mu', 0., 1.)
        rho = pm.HalfNormal('rho')

        # Priors are posteriors from previous iteration
        sigma1 = from_posterior('sigma1', trace['sigma1'])
        sigma2 = from_posterior('sigma2', trace['sigma2'])
        sigma3 = from_posterior('sigma3', trace['sigma3'])
        sigma4 = from_posterior('sigma4', trace['sigma4'])
    
        like1 = pm.Normal('like1', mu, sigma1, observed=measures[j][0])
        like2 = pm.Normal('like2', mu, sigma2, observed=measures[j][1])
        like3 = pm.Normal('like3', mu, sigma3, observed=measures[j][2])
        like4 = pm.Normal('like4', mu, sigma4, observed=measures[j][3])
        
        trace = pm.sample(1000)

Unfortunately the code inside the loop doesn’t work. The error is:

ValueError: Bad initial energy: inf. The model might be misspecified.

I tried with:

for RV in model.basic_RVs:
    print(RV.name, RV.logp(model.test_point))

But I have not found the variable that gives me .inf. Do you have any idea? Thank you.
I hope I do not seem too insistent, I’m trying to learn as much as I can.