introduction

J'ai essayé de visualiser la propriété d'agression de la distribution de probabilité représentée par la théorie de la limitation du pôle central avec Python. La nature atrophique est principalement extraite du livre "Statistics for Officially Certified Statistical Test Level 1 of the Japan Statistical Society".

Limitation du pôle central

Supposons que la distribution de probabilité $ X_i (i = 1, \ cdots, n) $ suit une même distribution indépendante de moyenne $ \ mu $ et de variance $ \ sigma ^ 2 $. Lorsque $ \ frac {\ sum_ {i = 1} ^ nX_i} {n} $ est $ n \ to \ infty $, la moyenne est $ \ mu $ et la variance est $ \ frac {\ sigma ^ 2} {n} $ Suivez une distribution normale.

Expérience

Comme c'est un gros problème, nous allons réaliser l'expérience avec une distribution déformée. $ X_i (i = 1, \ cdots, n) $ fonction de densité

f(x) = 11 x^{10}\ \ \ (0\leq x\leq1)

Et.

L'histogramme bleu est la valeur numérique de la fonction de densité de $ \ frac {\ sum_ {i = 1} ^ nX_i} {n} $, et la ligne orange continue est la fonction de densité de la distribution normale qui converge théoriquement. ..

n=1
n=2
Si $ n = 10 $

Cliquez ici pour le code source Python lorsque $ n = 1 $.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

#Moyenne de la fonction de densité de probabilité 11/12,Dispersion 11/13-(11/12)^2
def f(x):
    return 11 * (x ** 10)
mu = 11 / 12
sigma2 = 11 / 13 - (11 / 12) ** 2
#Fonction de distribution cumulative
def F(x):
    return x ** 11
#Fonction inverse de la fonction de distribution cumulative
def F_inv(x):
    return x ** (1 / 11)

n = 1
xmin = 0.6
xmax = 1.2
meanX = []
for _ in range(100000):
    meanX.append(np.sum(F_inv(np.random.rand(n))) / n)
plt.hist(meanX, bins=200, density=True)
s = np.linspace(xmin, xmax, 100)
plt.plot(s, norm.pdf(s, mu, (sigma2 / n) ** 0.5), linewidth=4)
plt.xlim(xmin, xmax)

Quantité de test statistique du test de conformité (1)

Supposons que $ (X_1, X_2, \ cdots, X_m) $ suit une distribution polynomiale avec $ n $ essais et $ (p_1, p_2, \ cdots, p_m) $ probabilité.

\sum_{i=1}^m\frac{(X_i-np_i)^2}{np_i}

Suit la distribution $ \ chi ^ 2 $ de $ (m-1) $ degrés de liberté lorsque $ n \ à \ infty $.

Expérience

Avec une distribution quaternaire avec une probabilité de $ \ big (\ frac {1} {16}, \ frac {1} {4}, \ frac {1} {4}, \ frac {7} {16} \ big) $ Pense. Obtenu numériquement

\sum_{i=1}^4\frac{(X_i-np_i)^2}{np_i}

Est représentée dans un histogramme bleu, et la distribution $ \ chi ^ 2 $ avec un degré de liberté 3 qui est théoriquement convergé est représentée par une ligne orange continue.

n=4
n=16
n=400

Cliquez ici pour le code source Python lorsque $ n = 4 $.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2
p = [1/16, 1/4, 1/4, 7/16]
n = 4
xmin = 0
xmax = 15
xx = []
for _ in range(100000):
    r = np.random.rand(1, n)
    x = [0] * 4
    x[0] = np.sum(r < sum(p[:1]))
    x[1] = np.sum(np.logical_and(sum(p[:1]) <= r, r < sum(p[:2])))
    x[2] = np.sum(np.logical_and(sum(p[:2]) <= r, r < sum(p[:3])))
    x[3] = np.sum(sum(p[:3]) <= r)
    xx.append(sum([(x[i] - n * p[i]) ** 2 / (n * p[i]) for i in range(4)]))
plt.hist(xx, bins=int(max(xx))*5, density=True)
s = np.linspace(xmin, xmax, 100)
plt.plot(s, chi2.pdf(s, 3), linewidth=4)
plt.xlim(xmin, xmax)

Quantité de test statistique du test de conformité (2)

La vraie probabilité $ p_i (i = 1,2, \ cdots, m) $ est inconnue, mais $ p_i $ est un paramètre de dimension $ l $ $ \ boldsymbol {\ theta} (l \ leq m-2) $ Supposons que vous sachiez que vous pouvez l'exprimer avec. Lorsque l'estimation la plus probable de $ p_i $ est $ \ hat {p} _i $

\sum_{i=1}^m\frac{(X_i-n\hat{p}_i)^2}{n\hat{p}_i}

Suit la distribution $ \ chi ^ 2 $ de $ (m-l-1) $ liberté quand $ n \ to \ infty $.

Expérience

Avec une distribution quaternaire avec une probabilité de $ \ big (\ frac {1} {16}, \ frac {1} {4}, \ frac {1} {4}, \ frac {7} {16} \ big) $ Pense. Le vrai $ p $ est inconnu, mais seul $ p_2 = p_3 $ est connu. À ce stade, il peut être exprimé comme $ p = [q, r, r, 1-2r-q] $. Trouvez $ q, r $ par la méthode d'estimation la plus probable.

\sum_{i=1}^m\frac{(X_i-n\hat{p}_i)^2}{n\hat{p}_i}

Est montré dans un histogramme bleu, et la distribution $ \ chi ^ 2 $ avec un degré de liberté 1 qui est théoriquement convergé est représentée en orange.

n=4
n=100
n=10000

Cliquez ici pour le code source Python lorsque $ n = 4 $.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2
n = 4
xmin = 0
xmax = 3
xx = []
for _ in range(100000):
    r = np.random.rand(1, n)
    x = [0] * 4
    x[0] = np.sum(r < sum(p[:1]))
    x[1] = np.sum(np.logical_and(sum(p[:1]) <= r, r < sum(p[:2])))
    x[2] = np.sum(np.logical_and(sum(p[:2]) <= r, r < sum(p[:3])))
    x[3] = np.sum(sum(p[:3]) <= r)
    p_ = [0] * 4
    p_[0] = x[0] / sum(x)
    p_[1] = (x[1] + x[2]) / (2 * sum(x))
    p_[2] = (x[1] + x[2]) / (2 * sum(x))
    p_[3] = x[3] / sum(x)  
    xx.append(sum([(x[i] - n * p_[i]) ** 2 / (n * p[i]) for i in range(4)]))
plt.hist(xx, bins=int(max(xx))*20, density=True)
s = np.linspace(xmin, xmax, 100)
plt.plot(s, chi2.pdf(s, 1), linewidth=4)
plt.xlim(xmin, xmax)

Normalité la plus proche de l'estimation la plus probable (partie 1)

Pour la variable stochastique $ X_i (i = 1, \ cdots, n) $ caractérisée par le paramètre $ \ theta $, $ \ theta_0 $ est le vrai paramètre, $ \ hat \ theta $ est l'estimation la plus probable, information de Fisher Montant $ J_n (\ theta) $

J_n(\theta)=E_\theta\Big[\Big(\frac{\delta}{\delta\theta}\log f(X_1,...,X_n;\theta)\Big)^2\Big]

Et. Lorsque $ \ hat \ theta $ vaut $ n \ to \ infty $, il suit une distribution normale de moyenne $ \ theta_0 $ et de variance $ J_n (\ theta_0) ^ {-1} $.

Expérience

$ X_i (i = 1, \ cdots, n) $ fonction de densité

f(x) = (\theta+1) x^{\theta}\ \ \ (0\leq x\leq1)

Et. Soit 10 le vrai $ \ theta $.

La distribution $ \ hat \ theta $ obtenue expérimentalement est représentée en bleu et la distribution normale théoriquement convergente est représentée en orange.

n=1
n=10
n=100

Cliquez ici pour le code source Python lorsque $ n = 1 $.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
#Moyenne de la fonction de densité de probabilité 11/12,Dispersion 11/13-(11/12)^2
def f(x):
    return 11 * (x ** 10)
mu = 11 / 12
sigma2 = 11 / 13 - (11 / 12) ** 2
#Fonction de distribution cumulative
def F(x):
    return x ** 11
#Fonction inverse de la fonction de distribution cumulative
def F_inv(x):
    return x ** (1 / 11)

n = 1
theta_min = -20
theta_max = 40
theta = []
for _ in range(100000):
    x = F_inv(np.random.rand(n))
    theta.append(- n / sum(np.log(x)) - 1)
theta = np.array(theta)

#histogramme thêta
th = np.linspace(theta_min, theta_max, 100)
hi = np.array([np.sum(np.logical_and(th[i] < theta, theta <= th[i + 1])) for i in range(100 - 1)]) / (len(theta) * (th[1] - th[0]))
plt.bar((th[:-1] + th[1:]) / 2, hi, width=(th[1] - th[0]))
#distribution normale
s = np.linspace(theta_min, theta_max, 100)
plt.plot(s, norm.pdf(s, 10, (11 ** 2 / n) ** 0.5), 'tab:orange', linewidth=4)
plt.xlim(theta_min, theta_max)
plt.ylim(0.0, 0.1)

Normalité la plus proche de l'estimation la plus probable (partie 2)

Soit $ g $ une fonction qui peut être différenciée par $ \ theta_0 $. Lorsque $ g (\ hat \ theta) $ est $ n \ to \ infty $, moyenne $ g (\ theta_0) $, variance $ g ^ \ prime (\ theta_0) ^ 2 J_n (\ theta_0) ^ {-1 } Suivez la distribution normale de $. Cependant, les définitions de $ \ theta $, $ \ hat \ theta $ et $ \ theta_0 $ sont les mêmes que dans la partie 1.

Expérience

Continuez l'expérience de la partie 1. Définissez $ g $ avec la formule suivante.

g(\theta)=\frac{1}{\theta}

La distribution $ g (\ hat \ theta) $ obtenue expérimentalement est représentée en bleu, et la distribution normale théoriquement convergente est représentée en orange.

n=1
n=10
n=100

Cliquez ici pour le code source Python lorsque $ n = 1 $.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

#Moyenne de la fonction de densité de probabilité 11/12,Dispersion 11/13-(11/12)^2
def f(x):
    return 11 * (x ** 10)
mu = 11 / 12
sigma2 = 11 / 13 - (11 / 12) ** 2
#Fonction de distribution cumulative
def F(x):
    return x ** 11
#Fonction inverse de la fonction de distribution cumulative
def F_inv(x):
    return x ** (1 / 11)

n = 1
theta_min = -0.2
theta_max = 0.6
theta = []
for _ in range(100000):
    x = F_inv(np.random.rand(n))
    theta.append(1 / (- n / sum(np.log(x)) - 1))
theta = np.array(theta)
th = np.linspace(theta_min, theta_max, 100)
hi = np.array([np.sum(np.logical_and(th[i] < theta, theta <= th[i + 1])) for i in range(100 - 1)]) / (len(theta) * (th[1] - th[0]))
plt.bar((th[:-1] + th[1:]) / 2, hi, width=(th[1] - th[0]))
s = np.linspace(theta_min, theta_max, 100)
plt.plot(s, norm.pdf(s, 1/10, (11 ** 2 / n / 10 ** 4) ** 0.5), 'tab:orange', linewidth=4)
plt.xlim(theta_min, theta_max)

Test statistique du test du rapport de vraisemblance

Supposons qu'il existe une distribution de probabilité $ f (x; \ theta) $ caractérisée par le paramètre $ \ theta $. Problème de test d'hypothèse

H_0:\theta\in\Theta_0 vs. $H_1:\theta\in\Theta_1 $

Dans, le rapport de vraisemblance $ L $ est défini par l'équation suivante.

L=\frac{\sup_{\theta\in\Theta_1} f(x_1,\cdots,x_n;\theta)}{\sup_{\theta\in\Theta_0} f(x_1,\cdots,x_n;\theta)}

Sous l'hypothèse nulle, $ 2 \ log L $ suit une distribution $ \ chi ^ 2 $ avec $ p $ liberté quand $ n \ to \ infty $. Cependant, $ p ＝ \ dim (\ Theta_1) - \ dim (\ Theta_0) $

Expérience

f(x;\theta) = (\theta+1) x^{\theta}\ \ \ (0\leq x\leq1)

En tant que problème de test d'hypothèse

H_0:\theta\=10 vs. $H_1\neq10 $

Et.

Sous l'hypothèse nulle, la fonction de densité de $ 2 \ log L $ est obtenue numériquement dans l'histogramme bleu, et la fonction de densité de la distribution $ \ chi ^ 2 $, qui est censée converger en théorie, est la ligne orange pleine.

Si $ n = 1 $
Si $ n = 10 $

Cliquez ici pour le code source lorsque $ n = 1 $.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2
#Distribution vraie: moyenne de la fonction de densité de probabilité 11/12,Dispersion 11/13-(11/12)^2
def f(x):
    return 11 * (x ** 10)

#Modèle paramétrique
def f(x, theta):
    return (theta + 1) * (x ** theta)

#Fonction de distribution cumulative
def F(x):
    return x ** 11

#Fonction inverse de la fonction de distribution cumulative
def F_inv(x):
    return x ** (1 / 11)

n = 1
l_min = 0
l_max = 5
l = []
for _ in range(100000):
    x = F_inv(np.random.rand(n))
    theta = - n / sum(np.log(x)) - 1  #Estimation la plus probable
    l.append(2 * np.log(np.prod(f(x, theta)) / np.prod(f(x, 10))))
l = np.array(l)
#l histogramme
th = np.linspace(l_min, l_max, 100)
hi = np.array([np.sum(np.logical_and(th[i] < l, l <= th[i + 1])) for i in range(100 - 1)]) / (len(l) * (th[1] - th[0]))
plt.bar((th[:-1] + th[1:]) / 2, hi, width=(th[1] - th[0]))
#Distribution du chi carré
s = np.linspace(l_min, l_max, 100)
plt.plot(s, chi2.pdf(s, 1), 'tab:orange', linewidth=4)
plt.xlim(l_min, l_max)
plt.ylim(0, 2)

Vérifiez la nature atrophique de la distribution de probabilité en Python

introduction

Limitation du pôle central

Expérience

Quantité de test statistique du test de conformité (1)

Expérience

Quantité de test statistique du test de conformité (2)

Expérience

Normalité la plus proche de l'estimation la plus probable (partie 1)

Expérience

Normalité la plus proche de l'estimation la plus probable (partie 2)

Expérience

Test statistique du test du rapport de vraisemblance

Expérience