[PYTHON] Test d'hypothèse statistique bayésienne

Test d'hypothèse statistique bayésienne

Dans cet article, je présenterai une méthode de test d'hypothèse avec des statistiques bayésiennes, et enfin je testerai réellement en utilisant python. Je suis un débutant en statistiques bayésiennes et en python, veuillez donc commenter s'il y a des lacunes dans le contenu.


1. Qu'est-ce que le test d'hypothèse?

$$ En gros, un test d'hypothèse est une hypothèse selon laquelle un paramètre inconnu $ \ theta $ appartient à un certain intervalle $ I_i $. Est-il correct de mettre en place une hypothèse $ H_i $ et d'adopter l'hypothèse $ H_i $, c'est-à-dire un paramètre inconnu? C'est une méthode pour analyser si l'on peut dire que $ \ theta $ appartient à l'intervalle $ I_i $.

Les trois méthodes suivantes sont typiques des statistiques bayésiennes.

Sous $$, la fonction de probabilité de la variable de probabilité $ X $ est représentée par $ p (X) $, et la probabilité que l'événement $ X $ soit représenté par $ P (X) $. De plus, en utilisant $ \ Omega $ comme espace d'échantillon, $ A \ bot B $ indique que $ A \ cap B = \ phi $ est valable pour l'ensemble $ A, B $.

2. Test par probabilité postérieure

L'idée du test en évaluant la probabilité ex post facto est très simple. Par exemple, considérons l'hypothèse suivante.

\left\{
\begin{split}
    H_0: &q \in I_0 \\
    H_1: &q \in I_0 \\
    &\vdots         \\
    H_k: &q \in I_k
\end{split}
\right.

$$ (Cependant, notez que $ I_i \ subset \ Omega $ et $ I_i $ peuvent partager de l'espace entre eux.)

$$ A ce moment, lorsque l'événement $ X $ est réalisé, les probabilités a posteriori que l'hypothèse $ H_i (i = 0, 1, \ dots k) $ tient sont différentes.

 P(H_i|X) = P(\theta \in I_i|X) = \int_{I_i} p(\theta|X) d\theta

Ce sera.

Si $ P (H_i | X) $ est assez grand (proche de 1), l'hypothèse $ H_i $ peut être adoptée, et si elle est petite (proche de 0), l'hypothèse $ H_i $ peut être rejetée.

3. Test par rapport de cotes postérieur

Vient ensuite un test utilisant le rapport de cotes postérieur. Vous pouvez faire autant d'hypothèses que vous le souhaitez dans le test de probabilité postérieure, mais seules deux hypothèses sont utilisées dans le test du rapport de cotes postérieur.

Considérez l'hypothèse suivante.

\left\{
\begin{split}
    H_0: q \in I_0 \\
    H_1: q \in I_1
\end{split}
\right.

$$ Ici, $ I_0 $ et $ I_1 $ sont définis de telle sorte que $ I_0 \ bot I_1, I_0 \ cup I_1 = \ Omega $ est valable. En faisant cela, l'un de $ H_0 $ et $ H_1 $ peut être l'hypothèse nulle et l'autre l'hypothèse alternative.

Le rapport de cotes ex post facto est

\begin{split}
Rapport de cotes post-hoc&= \frac{P(H_0|X)}{P(H_1|X)} \\
			 &= \frac{P(\theta \in I_0|X)}{P(\theta \in I_1|X)} \\
			 &= \frac{\int_{I_0}p(\theta|X)d\theta}{\int_{I_1}p(\theta|X)d\theta}
\end{split}

Il est défini par $$. Si le rapport de cotes postérieur est suffisamment supérieur à 1, $ H_0 $ sera adopté, et s'il est suffisamment petit, $ H_1 $ sera adopté.

$ Ce qui précède est la méthode de test utilisant le rapport de cotes postérieur, mais il y a quelques problèmes étant donné que la distribution postérieure est affectée par la distribution antérieure. Si la distribution antérieure est extrêmement favorable à l'une ou l'autre hypothèse, la distribution postérieure peut ne pas être correctement testée pour refléter les effets de la distribution antérieure. Par exemple, en réalitéH_1$Est correctH_0En supposant une distribution préalable extrêmement favorableP(H_0|X)Est plus grand qu'il ne l'est vraimentP(H_1|X)Est sous-estimé et le rapport de cotes postérieur peut être surestimé. Afin d'éviter cela, nous pouvons tester en utilisant le facteur bayésien.

4. Test factoriel bayésien

$$ Immédiatement, le facteur bayésien $ B_ {01} $ est défini par la formule suivante.

\begin{split} B_{01} &=Rapport de cotes post-hoc\rapport div pré-impair\\
          &= \frac{P(H_0|X)}{P(H_1|X)} \div \frac{P(H_0)}{P(H_1)} \\
          &=\frac{\int_{I_0}p(\theta|X)d\theta}{\int_{I_1}p(\theta|X)d\theta} \div \frac{\int_{I_0}p(\theta)d\theta}{\int_{I_1}p(\theta)d\theta}
\end{split}

Voyons maintenant pourquoi cette évaluation du facteur bayésien évite le problème du test de rapport de cotes postérieur.

Supposons que le réglage de la distribution antérieure $$ soit extrêmement avantageux à $ H_0 $, et que le rapport pré-cote est de 100 et le rapport post-cote est 20. À ce stade, $ H_0 $ sera adopté dans le test en fonction du rapport de cotes postérieur. Aussi, juste au cas où, si vous pensez à ce qui se passe avec le test par post-établissement, ce sera $ P (H_0 | X) \ simeq 0,95 $, et on peut dire que $ H_0 $ est également adopté ici. Cependant, si vous calculez le facteur bayésien $ B_ {01} $, vous obtenez $ B_ {01} = 0,2 $. La valeur du facteur Bayes est évaluée à l'aide des critères d'évaluation de Jeffreys indiqués ci-dessous. Puisqu'il est $ 1/10 \ le 0.2 \ le 1/3 $, dans ce cas, nous supportons pleinement l'hypothèse alternative $ H_1 $.

無題.png

[Source: https://www.slideshare.net/masarutokuoka/mcmc-35640699]

L'exemple ci-dessus est extrême, mais le test réel utilisant le facteur Bayes peut conduire à une conclusion différente du cas de l'utilisation du rapport de cotes postérieur, et le test utilisant le facteur Bayes est plus affecté par la distribution antérieure. On pense qu'il est difficile de tirer une fausse conclusion en réduisant la taille.

5. Testez en fait les exemples de données en utilisant python

J'ai en fait calculé la probabilité postérieure, le rapport de cotes postérieur et le facteur de Bayes et j'ai essayé le test d'hypothèse. Cette fois, par souci de simplicité, nous avons supposé une distribution bêta comme distribution a priori, et avons effectué un test d'hypothèse sur un modèle qui suit la distribution de Bernoulli.

Tout d'abord, définissez une fonction qui prend des exemples de données (données), des hyperparamètres de distribution antérieure (a0, b0) et un intervalle (I) comme arguments et renvoie la probabilité postérieure (post_prob), le rapport de cotes postérieur (post_odds_ratio) et le facteur bayésien (BF). Fait.

import scipy.stats as st

def ppb(data, a0, b0, I):
    n = data.size
    y = data.sum()
    pre_p0 = st.beta.cdf(I[1], a0, b0) - st.beta.cdf(I[0], a0, b0)
    pre_p1 = 1 - pre_p0
    post_p0 = st.beta.cdf(I[1], y + a0, n - y + b0) - st.beta.cdf(I[0], y + a0, n - y + b0)
    post_p1 = 1 - post_p0
    post_prob = post_p0
    post_odds_ratio = post_p0 / post_p1
    BF = post_odds_ratio * pre_p1 / pre_p0
    return post_prob, post_odds_ratio, BF

Nous avons ensuite généré aléatoirement des données à partir de la distribution de Bernoulli avec le paramètre p = 0,7 pour générer un échantillon d'une taille d'échantillon de 100.

import numpy as np

p = 0.7
n = 100
np.random.seed(0)
data = st.bernoulli.rvs(q, size=n)

Pour ces données, l'hyper paramètre

--A) a0 = b0 = 1 --B) a0 = 2, b0 = 5 --C) a0 = 3, b0 = 2

Nous avons analysé le cas de. L'intervalle (I) utilisé comme hypothèse nulle est I = [0,6, 0,8].

a0 = [1,2,3]
b0 = [1,5,2]
I = [0.6,0.8]
s = ['UNE','je','C']

for i in range(len(a0)):
    a, b, c = ppb(data, a0[i], b0[i], I)
    print('{}) \n probabilité postérieure: {:.5f} \n Odds ratio post-hoc: {:.5f} \n Facteur de Bayes: {:.5f}'.format(s[i],a,b,c))

Résultat d'exécution

UNE)
Probabilité post-hoc: 0.79632
Rapport de cotes post-hoc: 3.90973
Facteur Bayes: 15.63891
je)
Probabilité post-hoc: 0.93228
Rapport de cotes post-hoc: 13.76698
Facteur Bayes: 336.00387
C)
Probabilité post-hoc: 0.81888
Rapport de cotes post-hoc: 4.52127
Facteur Bayes: 8.62195

Ici, lorsque le test est effectué en référence aux critères d'évaluation du facteur bayésien par Jeffreys présentés précédemment, il est conclu que a) l'hypothèse nulle est fortement soutenue, b) l'hypothèse nulle est très fortement soutenue, et c) l'hypothèse nulle est suffisamment soutenue. A été fait. Puisque nous avons fixé la valeur réelle de p à 0,7, nous pouvons dire que nous avons tiré des conclusions raisonnables.

En outre, en regardant le graphique de distribution bêta ci-dessous, a) est le plus susceptible d'être affecté par la distribution précédente.

beta.png

L'image ci-dessous est un dessin réel du graphique de distribution postérieure. Certes, b) semble être plus influencé par la distribution antérieure que a) et c).

事後分布.png

De plus, lorsque seule la section d'hypothèse nulle a été modifiée en [0,65, 0,75], les résultats suivants ont été obtenus.

UNE)
Probabilité post-hoc: 0.34440
Rapport de cotes post-hoc: 0.52532
Facteur Bayes: 4.72784
je)
Probabilité post-hoc: 0.57232
Rapport de cotes post-hoc: 1.33818
Facteur Bayes: 74.33720
C)
Probabilité post-hoc: 0.36755
Rapport de cotes post-hoc: 0.58115
Facteur Bayes: 2.73401

Les conditions sont de plus en plus strictes et il est difficile de juger à partir de la probabilité postérieure et du rapport de cotes postérieur seuls, mais il semble que l'hypothèse nulle ne puisse être étayée dans les deux cas en utilisant le facteur bayésien.

Enfin, j'ai défini l'intervalle de l'hypothèse nulle comme [0,4, 0,6] et j'ai expérimenté quel type de conclusion peut être obtenu lorsque la valeur réelle n'est pas incluse dans l'intervalle.

UNE)
Probabilité post-hoc: 0.00019
Rapport de cotes post-hoc: 0.00019
Facteur Bayes: 0.00078
je)
Probabilité post-hoc: 0.00123
Rapport de cotes post-hoc: 0.00123
Facteur Bayes: 0.00515
C)
Probabilité post-hoc: 0.00020
Rapport de cotes post-hoc: 0.00020
Facteur Bayes: 0.00049
77

Il n'est pas nécessaire d'utiliser le facteur bayésien, mais dans les deux cas l'hypothèse nulle est correctement rejetée.

6. Enfin

J'ai supposé la distribution antérieure comme une distribution bêta et je n'ai traité que des modèles qui suivent la distribution de Bernoulli, mais j'aimerais traiter de modèles plus divers à l'avenir. Si vous avez des lacunes ou des erreurs, veuillez me le faire savoir. Merci d'avoir lu jusqu'au bout.


Les références

[Teruo Nakatsuma; Introduction aux statistiques bayésiennes avec Python, Asakura Shoten, 2019](https://www.amazon.co.jp/ by Python-Introduction to Bayesian Statistics-Practical Python Library-Nakazuma-Teruo / dp / 4254128983)

Rapport de recherche au MCMC (voir 2020-4-24)

Recommended Posts

Test d'hypothèse statistique bayésienne
Test d'hypothèse pour l'amélioration du produit
Test d'hypothèse et distribution de probabilité
tester