[PYTHON] Konzept des Bayes'schen Denkens (3) ... Berechnung der Änderungspunkte in der Anzahl der von PyMC3 empfangenen E-Mails
Aus Bayes'schen Überlegungen in Python
"In Python erlebte Bayes-Inferenz"
Hat sich die Anzahl der Nachrichten geändert?
Berechnen wir nun das Beispiel. Das "Python" -Skript ist in "PyMC2" im Buch geschrieben, aber nach meinem Verständnis basiert es auf dem "PyMC3" -Skript auf der offiziellen Website.
Hat sich die Anzahl der empfangenen Nachrichten geändert, wenn ein Benutzer jeden Tag eine Anzahl von Nachrichten empfängt? Das ist was es bedeutet. Holen Sie sich die Daten vom offiziellen "Github".
# -*- coding: utf-8 -*-
#Datenvisualisierung
from matplotlib import pyplot as plt
import numpy as np
count_data = np.loadtxt("data/txtdata.csv")
n_count_data = len(count_data)
plt.bar(np.arange(n_count_data), count_data, color="#348ABD")
plt.xlabel("Time (days)")
plt.ylabel("count of text-msgs received")
plt.title("Did the user's texting habits change over time?")
plt.xlim(0, n_count_data)
plt.show()
Hier ist die Grafik.
Sehen Sie das irgendwo als Veränderung? Sie können es nicht einfach anhand der Grafik erkennen. Was denken Sie
Konzept des Bayes'schen Denkens ... "Bayesismus und Wahrscheinlichkeitsverteilung"
Erstens ist dies die Anzahl der täglichen Empfänge, es handelt sich also um diskrete Wertdaten. Es hat also eine Poisson-Verteilung. Deshalb. Wenn die Anzahl der Nachrichten am $ i $ -Tag $ C_i $ beträgt,
C_i \sim \text{Poisson}(\lambda)
Was halten Sie von $ λ $ hier? Ich denke, dass sich $ λ $ in dieser Zeit irgendwo geändert hat. Beachten Sie, dass sich der Parameter $ λ $ am Tag $ τ $ geändert hat. Diese plötzliche Änderung wird als Schaltpunkt bezeichnet.
$ λ $ wird wie folgt ausgedrückt.
\lambda =
\begin{cases}
\lambda_1 & \text{if } t \lt \tau \cr
\lambda_2 & \text{if } t \ge \tau
\end{cases}
Wenn sich die Anzahl der empfangenen E-Mails nicht ändert, sollte dies $ \ lambda_1 = \ lambda_2 $ sein.
Um dieses $ λ $ zu betrachten, werden wir es mit einer Exponentialverteilung modellieren. (Weil $ λ $ ein stetiger Wert ist) Wenn der Parameter zu diesem Zeitpunkt $ α $ ist, wird er durch die folgende Formel ausgedrückt.
\begin{align}
&\lambda_1 \sim \text{Exp}( \alpha ) \\\
&\lambda_2 \sim \text{Exp}( \alpha )
\end{align}
Da der erwartete Wert der Exponentialverteilung die Umkehrung der Parameter ist, wird er wie folgt ausgedrückt.
\frac{1}{N}\sum_{i=0}^N \;C_i \approx E[\; \lambda \; |\; \alpha ] = \frac{1}{\alpha}
Dies schließt die Subjektivität in der vorherigen Verteilung nicht ein. Die Verwendung unterschiedlicher vorheriger Verteilungen für die beiden $ \ lambda_i $ drückt hier die Überzeugung aus, dass sich die Anzahl der Empfänge während des Beobachtungszeitraums irgendwo geändert hat.
Danach ist über $ τ $ zu denken, dass $ τ $ eine gleichmäßige Verteilung verwendet. Mit anderen Worten, der Glaube, dass jeder Tag der gleiche ist.
\begin{align}
& \tau \sim \text{DiscreteUniform(1,70) }\\\\
& \Rightarrow P( \tau = k ) = \frac{1}{70}
\end{align}
An dieser Stelle frage ich mich, wie man ein Programm mit python schreibt, aber danach werde ich es mit pymc3 schreiben.
import pymc3 as pm
import theano.tensor as tt
with pm.Model() as model:
alpha = 1.0/count_data.mean()
# Recall count_data is the
# variable that holds our txt counts
lambda_1 = pm.Exponential("lambda_1", alpha)
lambda_2 = pm.Exponential("lambda_2", alpha)
tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data - 1)
lambda_1 = pm.Exponential("lambda_1", alpha)
Erstellen Sie PyMC-Variablen, die $ \ lambda_1 $ und $ \ lambda_2 $ entsprechen.
tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data - 1)
Geben Sie $ \ tau $ mit einer gleichmäßigen Verteilung pm.DiscreteUniform an.
with model:
idx = np.arange(n_count_data) # Index
lambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)
with model:
observation = pm.Poisson("obs", lambda_, observed=count_data)
lambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)
Erstellen Sie nun ein neues lambda_.
Die Funktion switch () weist den Wert von lambda_1 oder lambda_2 als Wert von lambda_ zu, je nachdem auf welcher Seite von tau Sie sich befinden. Der Wert von "Lambda_" bis zu "Tau" ist "Lambda_1", und die Werte danach sind "Lambda_2".
observation = pm.Poisson("obs", lambda_, observed=count_data)
Dies bedeutet, dass die Daten "count_data" mit den von der Variablen "lambda_" bereitgestellten Variablen berechnet und beobachtet werden sollen.
with model:
step = pm.Metropolis()
trace = pm.sample(10000, tune=5000,step=step, cores=1)
lambda_1_samples = trace['lambda_1']
lambda_2_samples = trace['lambda_2']
tau_samples = trace['tau']
Die Details dieses Codes scheinen in Kapitel 3 des Buches zu sein. Dies ist der Lernschritt, der von einem Algorithmus namens Markov Chain Monte Carlo (MCMC) berechnet wird. Es gibt Tausende von stochastischen Variablen aus der posterioren Verteilung von $ \ lambda_1
Die Berechnung beginnt hier, aber da Multi-Core und Cuda nicht verwendet werden können (cores = 1), wird es lange dauern, aber die Antwort wird herauskommen.
Es ist eine grafische Darstellung des Ergebnisses.
Wie Sie den Ergebnissen entnehmen können
$ \ Tau $ beträgt oft 45 Tage und die Wahrscheinlichkeit beträgt 50%. Mit anderen Worten, ich kenne den Grund nicht, aber um den 45. Tag herum gab es einige Änderungen (einige Aktionen, die die Anzahl der Empfänge ändern ... aber dies ist nur der Person bekannt), und $ \ lambda $ hat sich geändert.
Schließlich werde ich den erwarteten Wert der Anzahl der Empfänge schreiben.
Auf diese Weise geben die Daten unter Verwendung der Bayes'schen Schätzung an, wann eine Änderung stattgefunden hat. Es ist jedoch zu beachten, dass der erwartete Wert der Anzahl der Empfänge tatsächlich eine Verteilung ist. Sie müssen auch selbst über die Ursache nachdenken.
Zum Schluss alle Skripte.
# -*- coding: utf-8 -*-
#Datenvisualisierung
from matplotlib import pyplot as plt
import numpy as np
count_data = np.loadtxt("data/txtdata.csv")
n_count_data = len(count_data)
plt.bar(np.arange(n_count_data), count_data, color="#348ABD")
plt.xlabel("Time (days)")
plt.ylabel("count of text-msgs received")
plt.title("Did the user's texting habits change over time?")
plt.xlim(0, n_count_data);
plt.show()
import pymc3 as pm
import theano.tensor as tt
with pm.Model() as model:
alpha = 1.0/count_data.mean()
# Recall count_data is the
# variable that holds our txt counts
lambda_1 = pm.Exponential("lambda_1", alpha)
lambda_2 = pm.Exponential("lambda_2", alpha)
tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data - 1)
with model:
idx = np.arange(n_count_data) # Index
lambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)
with model:
observation = pm.Poisson("obs", lambda_, observed=count_data)
### Mysterious code to be explained in Chapter 3.
with model:
step = pm.Metropolis()
trace = pm.sample(10000, tune=5000,step=step, cores=1)
lambda_1_samples = trace['lambda_1']
lambda_2_samples = trace['lambda_2']
tau_samples = trace['tau']
#histogram of the samples:
ax = plt.subplot(311)
ax.set_autoscaley_on(False)
plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85,
label="posterior of $\lambda_1$", color="#A60628", density=True)
plt.legend(loc="upper left")
plt.title(r"""Posterior distributions of the variables
$\lambda_1,\;\lambda_2,\;\tau$""")
plt.xlim([15, 30])
plt.xlabel("$\lambda_1$ value")
ax = plt.subplot(312)
ax.set_autoscaley_on(False)
plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.85,
label="posterior of $\lambda_2$", color="#7A68A6", density=True)
plt.legend(loc="upper left")
plt.xlim([15, 30])
plt.xlabel("$\lambda_2$ value")
plt.subplot(313)
w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples)
plt.hist(tau_samples, bins=n_count_data, alpha=1,
label=r"posterior of $\tau$",
color="#467821", weights=w, rwidth=2.)
plt.xticks(np.arange(n_count_data))
plt.legend(loc="upper left")
plt.ylim([0, .75])
plt.xlim([35, len(count_data)-20])
plt.xlabel(r"$\tau$ (in days)")
plt.ylabel("probability");
plt.show()
# tau_samples, lambda_1_samples, lambda_2_samples contain
# N samples from the corresponding posterior distribution
N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_count_data)
for day in range(0, n_count_data):
# ix is a bool index of all tau samples corresponding to
# the switchpoint occurring prior to value of 'day'
ix = day < tau_samples
# Each posterior sample corresponds to a value for tau.
# for each day, that value of tau indicates whether we're "before"
# (in the lambda1 "regime") or
# "after" (in the lambda2 "regime") the switchpoint.
# by taking the posterior sample of lambda1/2 accordingly, we can average
# over all samples to get an expected value for lambda on that day.
# As explained, the "message count" random variable is Poisson distributed,
# and therefore lambda (the poisson parameter) is the expected value of
# "message count".
expected_texts_per_day[day] = (lambda_1_samples[ix].sum()
+ lambda_2_samples[~ix].sum()) / N
plt.plot(range(n_count_data), expected_texts_per_day, lw=4, color="#E24A33",
label="expected number of text-messages received")
plt.xlim(0, n_count_data)
plt.xlabel("Day")
plt.ylabel("Expected # text-messages")
plt.title("Expected number of text-messages received")
plt.ylim(0, 60)
plt.bar(np.arange(len(count_data)), count_data, color="#348ABD", alpha=0.65,
label="observed texts per day")
plt.legend(loc="upper left")
plt.show()