Anpassungsmodell für zwei Normalverteilungen in PyMC

10

Da ich ein Softwareentwickler bin, der versucht, mehr Statistiken zu lernen, müssen Sie mir vergeben, bevor ich überhaupt anfange. Dies ist ein ernstes Neuland ...

Ich habe PyMC gelernt und einige wirklich (wirklich) einfache Beispiele durchgearbeitet . Ein Problem, bei dem ich nicht zur Arbeit kommen kann (und für das ich keine verwandten Beispiele finden kann), ist das Anpassen eines Modells an Daten, die aus zwei Normalverteilungen generiert wurden.

Angenommen, ich habe 1000 Werte. 500 generiert von a Normal(mean=100, stddev=20)und weitere 500 generiert von a Normal(mean=200, stddev=20).

Wenn ich ihnen ein Modell anpassen möchte, dh die beiden Mittelwerte und die einzelne Standardabweichung mithilfe von PyMC bestimmen. Ich weiß, es ist etwas in der Art von ...

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)

Das heißt, der Erzeugungsprozess ist Normal, aber mu ist einer von zwei Werten. Ich weiß nur nicht, wie ich die "Entscheidung" zwischen einem Wert von m1oder darstellen soll m2.

Vielleicht gehe ich einfach völlig falsch vor, um dies zu modellieren? Kann mich jemand auf ein Beispiel hinweisen? Ich kann BUGS und JAGS lesen, also ist wirklich alles in Ordnung.

mat kelcey
quelle

Antworten:

11

Sind Sie absolut sicher, dass die Hälfte von einer Distribution und die andere Hälfte von der anderen stammt? Wenn nicht, können wir den Anteil als Zufallsvariable modellieren (was sehr bayesianisch ist).

Folgendes würde ich tun, einige Tipps sind eingebettet.

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
quelle
2
Schamlose Werbung: Ich habe gerade einen Blog-Artikel über Bayes und pyMC geschrieben, buchstäblich 1 Minute bevor Sie dies gepostet haben, also lade ich Sie ein, es sich anzusehen . Die unglaubliche Kraft von Bayes - Teil 1
Cam.Davidson.Pilon
genial! Dieser Ansatz zur Vermischung der beiden Mittel ist genau das, was ich versucht habe, meinen Kopf herumzukriegen.
Mat Kelcey
Ich bin mir nicht sicher, ob ich den wahren Modellierungsvorteil vollständig verstehe, wenn ich sage, dass mean1 und mean2 normalverteilt statt einheitlich sind (das gleiche gilt wirklich für die Präzision, um ehrlich zu sein, ich verwende Gamma seit "jemand anderem"). Ich muss noch viel lernen :)
Mat Kelcey
Die Verwendung einer Uniform wie in Ihrem ursprünglichen Beispiel bedeutet, dass Sie mit absoluter Sicherheit wissen , dass der Mittelwert einen bestimmten Wert nicht überschreitet. Das ist etwas pathologisch. Es ist besser, eine Normale zu verwenden, da alle reellen Zahlen berücksichtigt werden können.
Cam.Davidson.Pilon
1
Die Wahl des Gammas hat einen mathematischen Grund. Das Gamma ist das konjugierte Prior der Präzision, siehe Tabelle hier
Cam.Davidson.Pilon
6

Einige Punkte im Zusammenhang mit der obigen Diskussion:

  1. Die Wahl zwischen diffusem Normal und Uniform ist ziemlich akademisch, es sei denn (a) Sie sind besorgt über die Konjugation. In diesem Fall würden Sie das Normal verwenden oder (b) es besteht eine vernünftige Wahrscheinlichkeit, dass der wahre Wert außerhalb der Endpunkte der Uniform liegt . Bei PyMC gibt es keinen Grund, sich über die Konjugation Gedanken zu machen, es sei denn, Sie möchten speziell einen Gibbs-Sampler verwenden.

  2. Ein Gamma ist eigentlich keine gute Wahl für einen Uninformativen vor einem Varianz- / Präzisionsparameter. Es kann informativer sein, als Sie denken. Eine bessere Wahl ist es, der Standardabweichung eine Uniform vorzuziehen und sie dann durch ein umgekehrtes Quadrat zu transformieren. Siehe Gelman 2006 für Details.

Fonnesbeck
quelle
1
ah fonnesbeck ist einer der kernentwickler von pymc! Können Sie uns ein Beispiel für die Codierung von Punkt 2 zeigen?
Cam.Davidson.Pilon
danke fonnesbeck und ja bitte! zu einem schnellen zB von Punkt 2 :)
mat kelcey
1
Ich vermute, Sie meinen etwas in der Art von ... gist.github.com/4404631 ?
Mat Kelcey
Ja genau. Sie können die Transformation etwas präziser durchführen:tau = std_dev**-2
Fonnesbeck
Was wäre der richtige Ort, um zu lesen, woher diese Beziehung zwischen Präzision und std_dev kommt?
user979