[PYTHON] Estimation bayésienne de chaque valeur attendue à partir d'une distribution bimodale

introduction

En utilisant la bibliothèque d'estimation bayésienne de Python PyMC3, de la distribution de probabilité discrète bimodale à la distribution de Poisson attendue $ \ lambda_1, \ lambda_2 $, et à chaque pic Estimez le rapport auquel il appartient.

Environnement d'exploitation et bibliothèque utilisés

Environnement d'exploitation

#OS
import platform
print(platform.platform())
#Windows-10-10.0.15063-SP0

#version python
print(platform.python_version())
#version 3.7.7

#Bibliothèque
import numpy as np #version 1.19.1
import pandas as pd #version 1.1.0
from scipy.stats import poisson #version 1.5.0
import matplotlib.pyplot as plt #version 3.2.2
import seaborn as sns #version 0.10.1
import pymc3 as pm #version 3.8
import theano.tensor as T #1.0.4

Essayer

En supposant que le nombre d'essais de Bernoulli est le même, considérons une distribution de probabilité discrète dans laquelle les comtés avec des nombres de succès attendus de $ \ lambda $ de 3 et 10 sont mélangés.

#Générez de manière aléatoire le nombre de succès lorsque la valeur attendue 3 est 900 et la valeur attendue 10 est essayée 1000 fois.
obs = np.hstack((poisson.rvs(3,size=900), poisson.rvs(10,size=1000)))

sns.countplot(obs)

Résultat de sortie (Je veux estimer le rapport avec la valeur attendue des deux pics de ce $ \ lambda $)

Considérons le modèle de PyMC3

Lorsque le comté avec la valeur attendue de 3 est A et le comté avec la valeur attendue de 10 est B, Puisqu'il n'est pas possible de savoir à l'avance à quel comté appartient la valeur observée $ n_x $, Supposons que la probabilité $ p (A) $ que la valeur observée $ n_x $ soit A suive une distribution uniforme de 0 à 1. Alors, la probabilité $ p (B) $ que la valeur observée $ n_x $ soit B est $ 1-p (A) $.

De plus, les valeurs attendues de A et B sont estimées en supposant que les véritables valeurs attendues de A et B sont également ** inconnues **. En supposant que la distribution précédente suit la distribution exponentielle, la valeur moyenne des valeurs observées est appliquée à $ \ lambda $ de la distribution exponentielle attendue $ \ frac {1} {\ lambda} $. En ce moment, je pense aussi qu'il y a deux pics. La distribution postérieure est estimée en utilisant la valeur obtenue à partir de la distribution antérieure comme valeur attendue mu et la valeur observée comme valeur générée aléatoirement plus tôt.

with pm.Model() as model:
    p1 = pm.Uniform('p1',0,1) #p(A)Probabilité de
    p2 = 1 - p1 #p(B)Probabilité de
    p = T.stack([p1, p2])
    
    #p1=0, p2=Préparez autant de variables catégorielles que 1 pour le nombre d'échantillons.
    assignment = pm.Categorical('assignment', p, shape=obs.shape[0], testval = np.random.randint(0, 2, obs.shape[0]))

    #Valeur attendue de la distribution antérieure= 1 /Défini comme une distribution exponentielle comme moyenne des valeurs observées
    lambda_ = pm.Exponential('lambda_', 1./obs.mean(), shape=2)
    lambda_i = pm.Deterministic('lambda_i', lambda_[assignment])
    #mu à lambda_Définissez la valeur générée sur observé
    observations = pm.Poisson('obs', mu=lambda_i, observed=obs)

    #chercher
    step1 = pm.Metropolis(vars=[p, lambda_]) #Valeurs discrètes définies sur Metropolis
    step2 = pm.ElemwiseCategorical(vars=[assignment])
    trace = pm.sample(5000, step=[step1, step2])

#Tracez les résultats
plt.figure(figsize=(12, 8))

plt.subplot(211)
plt.plot(trace['lambda_'][:,1], label='groupe A', c='r')
plt.plot(trace['lambda_'][:,0], label='Groupe B', c='c')
plt.ylim(2, 11)
plt.xlabel('Nombre de recherches')
plt.ylabel('Valeur attendue')
plt.title('Résultat estimé')
leg = plt.legend(loc='upper right')
leg.get_frame().set_alpha(0.7)

plt.subplot(212)
plt.plot(1 - trace['p1'], label='Pourcentage d'échantillons sélectionnés pour le comté A', c='r')
plt.legend(loc='upper right')
plt.xlabel('Nombre de recherches')
plt.ylabel('Pourcentage(%)')

plt.show()plt.figure(figsize=(12, 8))

#Affiche la valeur moyenne des recherches après la 1000e fois
trace_lambda_1 = trace['lambda_'][:,1][1000:].mean()
trace_lambda_2 = trace['lambda_'][:,0][1000:].mean()
p1 = (1 - trace['p1'][1000:]).mean()

print('Un résultat d'estimation de la valeur attendue du comté%.2f' % trace_lambda_1)
print('Résultat de l'estimation de la valeur attendue du groupe B%.2f' % trace_lambda_2)
print('Pourcentage du groupe A%.3f' % p1)

Un résultat de recherche de valeur attendue de comté 3.01
Résultat de recherche de valeur attendue du groupe B 9.84
Pourcentage du comté A 0.461

résultat

est devenu.

La distribution de Poisson attendue est de 3,01 et 9,84 La distribution de probabilité discrète dans laquelle l'échantillon consiste en un rapport de 46:54 est la suivante.

#0 dans la valeur attendue du résultat de l'estimation.001~0.Nombre de succès k prenant une probabilité entre 999
a_obs = np.arange(poisson.ppf(0.001, trace_lambda_1),
             poisson.ppf(0.999, trace_lambda_1))
b_obs = np.arange(poisson.ppf(0.001, trace_lambda_2),
             poisson.ppf(0.999, trace_lambda_2))
#Valeur attendue 3,0 à 10.001~0.Nombre de succès k prenant une probabilité entre 999
a_true = np.arange(poisson.ppf(0.001, 3),
             poisson.ppf(0.999, 3))
b_true = np.arange(poisson.ppf(0.001, 10),
             poisson.ppf(0.999, 10))
#Probabilité du nombre de succès k dans le résultat de l'estimation
A_obs = pd.DataFrame({'lambda': a_obs,
                 'percent': poisson.pmf(a_obs, trace_lambda_1) * (46 / 100)})
B_obs = pd.DataFrame({'lambda': b_obs,
                'percent': poisson.pmf(b_obs, trace_lambda_2) * (54 / 100)} )
#Valeur attendue 3,Probabilité de succès k à 10
A_true = pd.DataFrame({'lambda': a_true,
                 'percent': poisson.pmf(a_true, 3) * (9 / 19)})
B_true = pd.DataFrame({'lambda': b_true,
                'percent': poisson.pmf(b_true, 10) * (10 / 19)} )

#Comparez avec ceux générés aléatoirement
plt.figure(figsize=(20, 6) )
plt.subplot(121)
sns.barplot(data=pd.merge(A_obs,B_obs,how='outer', on=['lambda','percent']).groupby('lambda', as_index=False)['percent'].sum(),
            x='lambda', y='percent')
plt.title('Résultat estimé')
plt.ylabel('probabilité')
plt.xlabel('Nombre de succès k')
plt.xticks(np.int8(np.arange(0,21)),np.int8(np.arange(0,21)))
plt.subplot(122)
sns.barplot(data=pd.merge(A_true,B_true,how='outer', on=['lambda','percent']).groupby('lambda', as_index=False)['percent'].sum(),
            x='lambda', y='percent')
plt.title('Valeur attendue 3,Valeur attendue 10 ratio 9:10 distribution de probabilité')
plt.ylabel('probabilité')
plt.xlabel('Nombre de succès k')
plt.xticks(np.int8(np.arange(0,21)),np.int8(np.arange(0,21)))
plt.show()

Nous avons pu estimer la valeur attendue et le ratio de la distribution bimodale.