Ich habe einige Beobachtungsdaten, für die ich Parameter schätzen möchte, und ich dachte, es wäre eine gute Gelegenheit, PYMC3 auszuprobieren.

Meine Daten sind als eine Reihe von Datensätzen strukturiert. Jede Aufzeichnung enthält zwei Beobachtungen, die sich auf einen festgelegten Zeitraum von einer Stunde beziehen. Eine Beobachtung ist die Gesamtzahl der Ereignisse, die während der angegebenen Stunde auftreten. Die andere Beobachtung ist die Anzahl der Erfolge innerhalb dieses Zeitraums. So kann beispielsweise ein Datenpunkt angeben, dass in einem bestimmten Zeitraum von 1 Stunde insgesamt 1000 Ereignisse stattgefunden haben und dass von diesen 1000 100 Erfolge waren. In einem anderen Zeitraum kann es insgesamt 1000000 Ereignisse geben, von denen 120000 Erfolge sind. Die Varianz der Beobachtungen ist nicht konstant und hängt von der Gesamtzahl der Ereignisse ab. Zum Teil möchte ich diesen Effekt steuern und modellieren.

Mein erster Schritt dabei, um die zugrunde liegende Erfolgsquote abzuschätzen. Ich habe den folgenden Code vorbereitet, der die Situation nachahmen soll, indem zwei Sätze von "beobachteten" Daten bereitgestellt werden, indem sie mit scipy generiert werden. Es funktioniert jedoch nicht richtig.
Was ich davon erwarte, ist:

loss_lambda_factor ist ungefähr 0.1total_lambda (und total_lambda_mu) ist ungefähr 120.

Stattdessen konvergiert das Modell sehr schnell, aber auf eine unerwartete Antwort.

total_lambda und total_lambda_mu sind jeweils scharfe Peaks um 5e5.loss_lambda_factor ist ungefähr 0.

Das Traceplot (das ich aufgrund einer Reputation unter 10 nicht posten kann) ist ziemlich uninteressant - schnelle Konvergenz und scharfe Spitzen bei Zahlen, die nicht den Eingabedaten entsprechen. Ich bin gespannt, ob an meiner Herangehensweise etwas grundlegend falsch ist. Wie sollte der folgende Code geändert werden, um das richtige / erwartete Ergebnis zu erzielen?

from pymc import Model, Uniform, Normal, Poisson, Metropolis, traceplot 
from pymc import sample 
import scipy.stats

totalRates = scipy.stats.norm(loc=120, scale=20).rvs(size=10000)
totalCounts = scipy.stats.poisson.rvs(mu=totalRates) 
successRate = 0.1*totalRates 
successCounts = scipy.stats.poisson.rvs(mu=successRate) 

with Model() as success_model: 
    total_lambda_tau= Uniform('total_lambda_tau', lower=0, upper=100000)
    total_lambda_mu = Uniform('total_lambda_mu', lower=0, upper=1000000)
    total_lambda = Normal('total_lambda', mu=total_lambda_mu, tau=total_lambda_tau)
    total = Poisson('total', mu=total_lambda, observed=totalCounts) 

    loss_lambda_factor = Uniform('loss_lambda_factor', lower=0, upper=1)
    success_rate = Poisson('success_rate', mu=total_lambda*loss_lambda_factor, observed=successCounts) 

with success_model: 
    step =  Metropolis() 
    success_samples = sample(20000, step) #, start)


plt.figure(figsize=(10, 10)) 
_ = traceplot(success_samples)