Estimación del sesgo de la moneda a partir de múltiples observaciones del adversario

Question

Tengo un conjunto de monedas sesgadas $N$ donde cada moneda tiene un sesgo diferente $\theta_i$ (probabilidad de que haya cabezas para $i$ La moneda es $\theta_i$ ). Este sesgo es desconocido y debe ser estimado.

Hay $K$ expertos que, independientemente unos de otros, lanzan cada moneda exactamente una vez, observan el valor y lo comunican. Por desgracia, los expertos no son fieles, con la probabilidad $p_j$ El $j$ El experto informará de la cara por la cruz, o de la cruz por la cara. Se sabe que la probabilidad de que salga cara es inferior a 0,5 para todos los expertos y es necesario estimarla.

¿Cuál sería el modelo adecuado para este problema? Soy nuevo en la programación probabilística y cualquier consejo sobre la lectura, la experimentación sería genial.

Estoy familiarizado con PyMC3, algo de JAGS y muy poco con Stan.

Answer 1

He ideado un modelo que estima el sesgo de la moneda y las probabilidades de lanzamiento de tres adversarios.

#
# Estimate coin bias given multiple observations from adversarial experts.
#

import sys

import pymc3 as pm
import numpy as np
import theano.tensor as tt
import matplotlib.pyplot as plt

from scipy.optimize import fmin_powell
from scipy.stats.distributions import bernoulli

p_true = 0.1                                # true coin bias (probability of heads)
a_true = np.array( [ 0.1, 0.2, 0.3 ] )      # noise (flipping probability)

N = 1000        # number of coin observations
K = 3           # number of experts

# generate reference data
X = bernoulli.rvs( size=N, p=p_true )

# corrupt data with noise
Y = np.zeros( (K,N) )
for k in range( K ):
    Y[k,:] = X
    flip_or_not = bernoulli.rvs( size=N, p=a_true[k] )
    for i in range( N ):
        if flip_or_not[i] == 1:
            if Y[k,i] == 1:
                Y[k,i] = 0
            else:
                Y[k,i] = 1

model = pm.Model()

with model:
    alpha0 = pm.HalfCauchy( 'alpha0', beta=1 )
    beta0 = pm.HalfCauchy( 'beta0', beta=1 )

    p = pm.Beta( 'p', alpha=alpha0, beta=beta0 )
    a = pm.Uniform( 'a', lower=0, upper=0.5, shape=K )

    q = a + p - 2 * a * p

    y_hidden = pm.Bernoulli('y_hidden_' + str(k), p=p, shape=N)

    for k in range( K ):
        y_obs = pm.Bernoulli( 'y_obs_' + str(k), p=q[k], observed=Y[k,:] )
        pot = pm.Potential( 'pot_' + str(k), -1000000 * ( tt.sum( (y_hidden - y_obs)**2 ) ) )

    #for i in range( N ):
    #    potential = pm.Potential( 'potential_' + str(i), -1000*(y_hidden[i] - y_obs[i]) ** 2 )

    start = pm.find_MAP()
    step = pm.NUTS( scaling=start )
    trace = pm.sample( 1000, start=start, step=step )

varnames = [ 'a', 'p' ]

pm.traceplot( trace, varnames=varnames )
#pm.autocorrplot( trace )
plt.show()