[PYTHON] Analyse par raisonnement bayésien (1) ... Quel est le meilleur, A ou B?
Exemple d'analyse du raisonnement bayésien;: Quel est le meilleur, A ou B?
Du raisonnement bayésien expérimenté en Python
"Inférence de Bayes expérimentée en Python"
Quel est le meilleur site A ou site B dans le raisonnement bayésien? C'est le problème.
Question: Quelle est la vraie conversion pour le site A?
Soit $ p_A $ la probabilité qu'un utilisateur qui voit le site A mène finalement à une action rentable (c'est ce qu'on appelle la conversion) telle que la demande de matériel ou l'achat. Supposons que $ N $ personnes regardent le site A et que $ n $ personnes mènent à des conversions. Puis
p_A=n/N
Vous pourriez penser, mais ce n'est pas le cas. Autrement dit, nous ne savons pas si $ p_A $ et $ n / N $ sont égaux. Il existe une différence entre la fréquence observée et la fréquence réelle, et la fréquence réelle peut être interprétée comme la probabilité qu'un événement se produise. En d'autres termes, depuis que j'ai vu le site A, je ne sais pas si cela a conduit à l'action.
Par exemple, en prenant un dé comme exemple, la fréquence réelle de lancer un dé et d'obtenir un 1 est de 1/6
Tout d'abord, nous devons considérer la distribution antérieure pour le nombre inconnu $ p_A $. Que pensez-vous de $ p_A $ lorsque vous n'avez pas les données? Je ne suis pas encore sûr, alors supposons que $ p_A $ suit une distribution uniforme de [0,1].
Donnez à la variable p une distribution uniforme de 0 à 1 avec @ pm.Uniform.
import pymc3 as pm
# The parameters are the bounds of the Uniform.
with pm.Model() as model:
p = pm.Uniform('p', lower=0, upper=1)
Par exemple, $ p_A = 0,05 $ et $ N = 1500 $ Simulons le nombre d'utilisateurs convertis (ont acheté ou demandé des matériaux en regardant le site) lorsque l'utilisateur a regardé le site A. En fait, je ne connais pas $ p_A $, mais je suppose que je le connais pour le moment et que je le simule pour créer des données d'observation.
Puisqu'il y a N personnes, [distribution Bernouy](https://ja.wikipedia.org/wiki/%E3%83%99%E3%83%AB%E3%83%] pour simuler N essais 8C% E3% 83% BC% E3% 82% A4% E5% 88% 86% E5% B8% 83) est utilisé. Puisque la distribution de Bernoulli est une distribution de probabilité pour les variables binaires (variables qui prennent 0 ou 1), c'est une bonne distribution de probabilité pour calculer s'il faut convertir ou non.
X\sim \text{Ber}(p)
À ce moment, $ X $ est de 1 $ avec une probabilité de $ p $ et de 0 $ avec une probabilité de 1 $-p $.
Le script python qui simule cela est
import scipy.stats as stats
import numpy as np
#set constants
p_true = 0.05 # remember, this is unknown.
N = 1500
# sample N Bernoulli random variables from Ber(0.05).
# each random variable has a 0.05 chance of being a 1.
# this is the data-generation step
occurrences = stats.bernoulli.rvs(p_true, size=N)
print(occurrences) # Remember: Python treats True == 1, and False == 0
print(np.sum(occurrences))
#[0 0 0 ... 0 0 0]
#67
# Occurrences.mean is equal to n/N.
print("What is the observed frequency in Group A? %.4f" % np.mean(occurrences))
print("Does this equal the true frequency? %s" % (np.mean(occurrences) == p_true))
#What is the observed frequency in Group A? 0.0447
#Does this equal the true frequency? False
Maintenant que nous avons les données observées (occurrences), nous définissons les données observées dans la variable obs de PyMC et exécutons l'algorithme d'inférence pour créer un graphique de la distribution postérieure de $ p_A $ inconnu.
#include the observations, which are Bernoulli
with model:
obs = pm.Bernoulli("obs", p, observed=occurrences)
# To be explained in chapter 3
step = pm.Metropolis()
trace = pm.sample(18000, step=step,cores =1)
burned_trace = trace[1000:]
import matplotlib.pyplot as plt
plt.title("Posterior distribution of $p_A$, the true effectiveness of site A")
plt.vlines(p_true, 0, 90, linestyle="--",color="r", label="true $p_A$ (unknown)")
plt.hist(burned_trace["p"], bins=25, histtype="bar",density=True)
plt.legend()
plt.show()
Le graphique résultant ressemble à ceci.
Cette fois, le graphique de la post-probabilité est comme ceci, tandis que la pré-probabilité de la valeur vraie est de 0,5. Si vous répétez ceci, vous obtiendrez un graphique légèrement différent à chaque fois.
Question: En comparant le site A et le site B, lequel est le meilleur?
Si le site B est calculé de la même manière, la probabilité postérieure de $ p_B $ peut être obtenue. J'ai donc décidé de calculer $ delta = p_A-p_B $.
Comme $ p_A $, $ p_B $ ne connaît pas la vraie valeur. Supposons que $ p_B = 0,4 $. Alors $ delta = 0,1 $. Calculez $ p_B $ et $ delta $ de la même manière que $ p_A $. J'ai tout calculé à nouveau.
# -*- coding: utf-8 -*-
import pymc3 as pm
import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt
#these two quantities are unknown to us.
true_p_A = 0.05
true_p_B = 0.04
#notice the unequal sample sizes -- no problem in Bayesian analysis.
N_A = 1500
N_B = 750
#generate some observations
observations_A = stats.bernoulli.rvs(true_p_A, size=N_A)
observations_B = stats.bernoulli.rvs(true_p_B, size=N_B)
print("Obs from Site A: ", observations_A[:30], "...")
print("Obs from Site B: ", observations_B[:30], "...")
#Obs from Site A: [0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] ...
#Obs from Site B: [0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] ...
print(np.mean(observations_A))
#0.04666666666666667
print(np.mean(observations_B))
#0.03866666666666667
# Set up the pymc3 model. Again assume Uniform priors for p_A and p_B.
with pm.Model() as model:
p_A = pm.Uniform("p_A", 0, 1)
p_B = pm.Uniform("p_B", 0, 1)
# Define the deterministic delta function. This is our unknown of interest.
delta = pm.Deterministic("delta", p_A - p_B)
# Set of observations, in this case we have two observation datasets.
obs_A = pm.Bernoulli("obs_A", p_A, observed=observations_A)
obs_B = pm.Bernoulli("obs_B", p_B, observed=observations_B)
# To be explained in chapter 3.
step = pm.Metropolis()
trace = pm.sample(20000, step=step,cores=1)
burned_trace=trace[1000:]
p_A_samples = burned_trace["p_A"]
p_B_samples = burned_trace["p_B"]
delta_samples = burned_trace["delta"]
#histogram of posteriors
ax = plt.subplot(311)
plt.xlim(0, .1)
plt.hist(p_A_samples, histtype='stepfilled', bins=25, alpha=0.85,
label="posterior of $p_A$", color="#A60628", density=True)
plt.vlines(true_p_A, 0, 80, linestyle="--", label="true $p_A$ (unknown)")
plt.legend(loc="upper right")
plt.title("Posterior distributions of $p_A$, $p_B$, and delta unknowns")
ax = plt.subplot(312)
plt.xlim(0, .1)
plt.hist(p_B_samples, histtype='stepfilled', bins=25, alpha=0.85,
label="posterior of $p_B$", color="#467821", density=True)
plt.vlines(true_p_B, 0, 80, linestyle="--", label="true $p_B$ (unknown)")
plt.legend(loc="upper right")
ax = plt.subplot(313)
plt.hist(delta_samples, histtype='stepfilled', bins=30, alpha=0.85,
label="posterior of delta", color="#7A68A6", density=True)
plt.vlines(true_p_A - true_p_B, 0, 60, linestyle="--",
label="true delta (unknown)")
plt.vlines(0, 0, 60, color="black", alpha=0.2)
plt.legend(loc="upper right");
plt.show()
# Count the number of samples less than 0, i.e. the area under the curve
# before 0, represent the probability that site A is worse than site B.
print("Probability site A is WORSE than site B: %.3f" % \
np.mean(delta_samples < 0))
print("Probability site A is BETTER than site B: %.3f" % \
np.mean(delta_samples > 0))
#Probability site A is WORSE than site B: 0.201
#Probability site A is BETTER than site B: 0.799
Le graphique résultant ressemble à ceci.
On voit que la base de la distribution postérieure de $ p_B $ est légèrement plus large. C'est parce que le site B a moins de données. Cela indique que nous ne sommes pas sûrs de la valeur de $ p_B $ en raison du manque de données sur le site B.
Le troisième graphique, $ delta = p_A-p_B $, montre une valeur supérieure à 0. Cela signifie que le site A est plus réactif aux utilisateurs en tant que conversion.
Vous pouvez obtenir une meilleure compréhension en changeant les valeurs de $ p_A $, $ p_B $ et le nombre de données.
Recommended Posts