[PYTHON] Concept de raisonnement bayésien (3) ... Calcul des points de changement dans le nombre d'emails reçus par PyMC3

Du raisonnement bayésien expérimenté en Python

"Inférence de Bayes expérimentée en Python"

Le nombre de messages a-t-il changé?

Calculons maintenant l'exemple. Le script Python est écrit en PyMC2 dans le livre, mais dans le cadre de ma compréhension, il est écrit sur la base du script PyMC3 sur le site officiel.

Lorsqu'un utilisateur a reçu un certain nombre de messages chaque jour, ce nombre de messages reçus a-t-il changé? C'est ce que cela signifie. Récupérez les données du github officiel.

# -*- coding: utf-8 -*-
#Visualisation de données
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()

Voici le graphique.

Voyez-vous cela comme un changement quelque part? Vous ne pouvez pas le dire simplement en regardant le graphique. Qu'est-ce que tu penses?

Concept de raisonnement bayésien ... "Bayésisme et distribution de probabilité"

Premièrement, il s'agit du nombre de réceptions quotidiennes, il s'agit donc de données à valeur discrète. Il a donc une distribution de Poisson. Donc. Si le nombre de messages le jour $ i $ est $ C_i $,

C_i \sim \text{Poisson}(\lambda) 

Que pensez-vous de $ λ $ ici? Je pense que $ λ $ a peut-être changé quelque part pendant cette période. Considérez ceci que le paramètre $ λ $ a changé le jour $ τ $. Ce changement soudain est appelé un point de commutation.

$ λ $ est exprimé comme suit.

\lambda = 
\begin{cases}
\lambda_1  & \text{if } t \lt \tau \cr
\lambda_2 & \text{if } t \ge \tau
\end{cases}

Si le nombre d'e-mails reçus ne change pas, il devrait être $ \ lambda_1 = \ lambda_2 $.

Maintenant, pour considérer ce $ λ $, nous allons le modéliser avec une distribution exponentielle. (Parce que $ λ $ est une valeur continue) Si le paramètre à ce moment est $ α $, il est exprimé par la formule suivante.

\begin{align}
&\lambda_1 \sim \text{Exp}( \alpha ) \\\
&\lambda_2 \sim \text{Exp}( \alpha )
\end{align}

Puisque la valeur attendue de la distribution exponentielle est l'inverse des paramètres, elle est exprimée comme suit.

\frac{1}{N}\sum_{i=0}^N \;C_i \approx E[\; \lambda \; |\; \alpha ] = \frac{1}{\alpha} 

Cela n'inclut pas la subjectivité dans la distribution antérieure. Ici, l'utilisation de distributions antérieures différentes pour les deux $ \ lambda_i $ exprime la croyance que le nombre de réceptions a changé quelque part pendant la période d'observation.

Après cela, ce qu'il faut penser de $ τ $ est que $ τ $ utilise une distribution uniforme. En d'autres termes, la croyance que chaque jour est le même.

\begin{align}
& \tau \sim \text{DiscreteUniform(1,70) }\\\\
& \Rightarrow P( \tau = k ) = \frac{1}{70}
\end{align}

À ce stade, je me demande comment écrire un programme avec python, mais après cela, je l'écrirai avec pymc3.

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) Créez des variables PyMC correspondant à $ \ lambda_1 $ et $ \ lambda_2 $.

tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data - 1) Donnez $ \ tau $ avec une distribution uniforme pm.DiscreteUniform.

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) Créez maintenant un nouveau lambda_. La fonction switch () affecte la valeur de lambda_1 ou lambda_2 comme valeur de lambda_, selon le côté de tau sur lequel vous vous trouvez. La valeur de lambda_ jusqu'à tau est lambda_1, et les valeurs suivantes sont lambda_2.

observation = pm.Poisson("obs", lambda_, observed=count_data)

Cela signifie que les données «count_data» doivent être calculées et observées avec les variables fournies par la variable «lambda_».

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']

Les détails de ce code semblent être dans le chapitre 3 du livre. Il s'agit de l'étape d'apprentissage, qui est calculée par un algorithme appelé Markov Chain Monte Carlo (MCMC). Il renvoie des milliers de variables stochastiques de la distribution postérieure de $ \ lambda_1 $ \ lambda_2 $ \ tau $ appelées traces. En traçant cette variable stochastique, vous pouvez voir la forme de la distribution postérieure.

Le calcul commence ici, mais comme le multi-core et Cuda ne peuvent pas être utilisés (cores = 1), cela prendra du temps, mais la réponse sortira.

C'est un graphique du résultat.

Comme vous pouvez le voir dans les résultats \lambda_1=18 \lambda_2=23 \tau=45

$ \ Tau $ est souvent de 45 jours, et la probabilité est de 50%. En d'autres termes, je ne connais pas la raison, mais vers le 45e jour, il y a eu un changement (une action qui change le nombre de réceptions ... mais cela n'est connu que de la personne), et $ \ lambda $ a changé.

Enfin, j'écrirai la valeur attendue du nombre de réceptions.

De cette façon, en utilisant l'estimation bayésienne, les données nous indiquent quand il y a eu un changement. Cependant, il convient de noter que la valeur attendue du nombre de réceptions est en fait une distribution. Vous devez également penser vous-même à la cause.

Enfin, tous les scripts.

# -*- coding: utf-8 -*-
#Visualisation de données
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()