Modèle d'ajustement pour deux distributions normales dans PyMC

10

Étant donné que je suis un ingénieur logiciel essayant d'apprendre plus de statistiques, vous devrez me pardonner avant même de commencer, c'est un nouveau territoire sérieux ...

J'ai appris PyMC et travaillé à travers des exemples vraiment (vraiment) simples. Un problème pour lequel je ne peux pas travailler (et je ne trouve aucun exemple connexe) est l'ajustement d'un modèle aux données générées à partir de deux distributions normales.

Disons que j'ai 1000 valeurs; 500 générés à partir d'un Normal(mean=100, stddev=20)et 500 autres générés à partir d'un Normal(mean=200, stddev=20).

Si je veux leur adapter un modèle, c'est-à-dire déterminer les deux moyennes et l'écart-type unique, en utilisant PyMC. Je sais que c'est quelque chose comme ...

mean1 = Uniform('mean1', lower=0.0, upper=200.0)
mean2 = Uniform('mean2', lower=0.0, upper=200.0)
precision = Gamma('precision', alpha=0.1, beta=0.1)

data = read_data_from_file_or_whatever()

@deterministic(plot=False)
def mean(m1=mean1, m2=mean2):
    # but what goes here?

process = Normal('process', mu=mean, tau=precision, value=data, observed=True)

c'est-à-dire que le processus de génération est Normal, mais mu est l'une des deux valeurs. Je ne sais tout simplement pas comment représenter la "décision" entre le fait qu'une valeur provienne de m1ou m2.

Peut-être que je prends complètement la mauvaise approche pour modéliser cela? Quelqu'un peut-il me donner un exemple? Je peux lire les BUGS et JAGS donc tout va bien vraiment.

mat kelcey
la source

Réponses:

11

Êtes-vous absolument certain que la moitié provenait d'une distribution et l'autre de l'autre? Sinon, nous pouvons modéliser la proportion comme une variable aléatoire (ce qui est une chose très bayésienne à faire).

Voici ce que je ferais, quelques conseils sont intégrés.

from pymc import *

size = 10
p = Uniform( "p", 0 , 1) #this is the fraction that come from mean1 vs mean2

ber = Bernoulli( "ber", p = p, size = size) # produces 1 with proportion p.

precision = Gamma('precision', alpha=0.1, beta=0.1)

mean1 = Normal( "mean1", 0, 0.001 ) #better to use normals versus Uniforms (unless you are certain the value is  truncated at 0 and 200 
mean2 = Normal( "mean2", 0, 0.001 )

@deterministic
def mean( ber = ber, mean1 = mean1, mean2 = mean2):
    return ber*mean1 + (1-ber)*mean2


#generate some artificial data   
v = np.random.randint( 0, 2, size)
data = v*(10+ np.random.randn(size) ) + (1-v)*(-10 + np.random.randn(size ) )


obs = Normal( "obs", mean, precision, value = data, observed = True)

model = Model( {"p":p, "precision": precision, "mean1": mean1, "mean2":mean2, "obs":obs} )
Cam.Davidson.Pilon
la source
2
Promotion sans vergogne: je viens d'écrire un article de blog sur Bayes et pyMC littéralement 1 minute avant de poster ceci, donc je vous invite à le vérifier The Awesome Power of Bayes - Part 1
Cam.Davidson.Pilon
impressionnant! cette approche du mélange des deux moyens est exactement ce que j'essayais de comprendre.
mat kelcey
Je ne suis pas sûr de bien comprendre le véritable avantage de la modélisation de dire que la moyenne1 et la moyenne2 sont normalement distribuées au lieu d'uniformes (il en va vraiment de la précision pour être honnête, j'utilise Gamma depuis "quelqu'un d'autre"). J'ai beaucoup à apprendre :)
mat kelcey
L'utilisation d'un uniforme, comme dans votre exemple d'origine, implique que vous savez avec une certitude absolue que la moyenne ne dépasse pas une certaine valeur. C'est quelque peu pathologique. Il est préférable d'utiliser une normale, car elle permet de prendre en compte tous les nombres réels.
Cam.Davidson.Pilon
1
Le choix du gamma a une raison mathématique. Le gamma est le conjugué avant la précision, voir le tableau ici
Cam.Davidson.Pilon
6

Quelques points, liés à la discussion ci-dessus:

  1. Le choix de la normale diffuse par rapport à l'uniforme est assez académique, sauf si (a) vous vous inquiétez de la conjugaison, auquel cas vous utiliseriez la normale ou (b) il y a une chance raisonnable que la vraie valeur puisse être en dehors des points limites de l'uniforme . Avec PyMC, il n'y a aucune raison de s'inquiéter de la conjugaison, sauf si vous souhaitez spécifiquement utiliser un échantillonneur Gibbs.

  2. Un gamma n'est en fait pas un excellent choix pour un non informatif avant un paramètre de variance / précision. Cela peut finir par être plus informatif que vous ne le pensez. Un meilleur choix est de mettre un a priori uniforme sur l'écart-type, puis de le transformer par un carré inverse. Voir Gelman 2006 pour plus de détails.

fonnesbeck
la source
1
ah fonnesbeck est l'un des principaux développeurs de pymc! Pouvez-vous nous montrer un exemple de codage du point 2?
Cam.Davidson.Pilon
merci fonnesbeck et, oui s'il vous plait! à un exemple rapide du point 2 :)
mat kelcey
1
en fait, je suppose que vous voulez dire quelque chose du genre ... gist.github.com/4404631 ?
mat kelcey
Oui, exactement. Vous pouvez faire la transformation un peu plus concise:tau = std_dev**-2
fonnesbeck
quel serait le bon endroit pour lire d'où vient cette relation entre la précision et std_dev?
user979