Überprüfen Sie die atrophische Natur der Wahrscheinlichkeitsverteilung in Python

Einführung

Ich habe versucht, die Aggressionseigenschaft der Wahrscheinlichkeitsverteilung zu visualisieren, die durch die zentrale Polbegrenzungstheorie mit Python dargestellt wird. Die atrophische Natur wird hauptsächlich aus dem Buch "Statistik für offiziell zertifizierte statistische Teststufe 1 der Japan Statistical Society" entnommen.

Zentralpolbegrenzung

Angenommen, die Wahrscheinlichkeitsverteilung $ X_i (i = 1, \ cdots, n) $ folgt einer unabhängigen gleichen Verteilung mit dem Mittelwert $ \ mu $ und der Varianz $ \ sigma ^ 2 $. Wenn $ \ frac {\ sum_ {i = 1} ^ nX_i} {n} $ $ n \ bis \ infty $ ist, ist der Durchschnitt $ \ mu $ und die Varianz ist $ \ frac {\ sigma ^ 2} {n} $ Folgen Sie einer Normalverteilung.

Experiment

Da es eine große Sache ist, werden wir das Experiment mit einer verzerrten Verteilung durchführen. $ X_i (i = 1, \ cdots, n) $ Dichtefunktion

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

Und.

Das blaue Histogramm ist der numerische Wert der Dichtefunktion von $ \ frac {\ sum_ {i = 1} ^ nX_i} {n} $, und die durchgezogene orange Linie ist die Dichtefunktion der Normalverteilung, die theoretisch konvergiert. ..

• n=1

• n=2

• Wenn $ n = 10 $

Klicken Sie hier für Python-Quellcode, wenn $ n = 1 $.

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

#Durchschnitt der Wahrscheinlichkeitsdichtefunktion 11/12,Dispersion 11/13-(11/12)^2
def f(x):
    return 11 * (x ** 10)
mu = 11 / 12
sigma2 = 11 / 13 - (11 / 12) ** 2
#Verteilungsfunktion
def F(x):
    return x ** 11
#Inverse Funktion der kumulativen Verteilungsfunktion
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)

Statistische Testmenge der Konformitätsprüfung (1)

Angenommen, $ (X_1, X_2, \ cdots, X_m) $ folgt einer Polynomverteilung mit $ n $ Versuchen und $ (p_1, p_2, \ cdots, p_m) $ Wahrscheinlichkeit.

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

Folgt der $ \ chi ^ 2 $ -Verteilung von $ (m-1) $ Freiheitsgraden, wenn $ n \ bis \ infty $.

Experiment

Mit einer quaternären Verteilung von $ \ big (\ frac {1} {16}, \ frac {1} {4}, \ frac {1} {4}, \ frac {7} {16} \ big) Überlegen. Numerisch erhalten

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

Wird im blauen Histogramm angezeigt, und die theoretisch konvergierte Verteilung von $ \ chi ^ 2 $ mit 3 Freiheitsgraden wird durch die durchgezogene orange Linie angezeigt.

• n=4

• n=16

• n=400

Klicken Sie hier für den Python-Quellcode, wenn $ 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)

Statistische Testmenge der Konformitätsprüfung (2)

Die wahre Wahrscheinlichkeit $ p_i (i = 1,2, \ cdots, m) $ ist unbekannt, aber $ p_i $ ist ein $ l $ -Dimensionsparameter $ \ boldsymbol {\ theta} (l \ leq m-2) $ Angenommen, Sie wissen, dass Sie es mit ausdrücken können. Wenn die wahrscheinlichste Schätzung von $ p_i $ $ \ hat {p} _i $ ist

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

Folgt der $ \ chi ^ 2 $ -Verteilung von $ (m-l-1) $ Freiheitsgraden, wenn $ n \ bis \ infty $.

Experiment

Mit einer quaternären Verteilung von $ \ big (\ frac {1} {16}, \ frac {1} {4}, \ frac {1} {4}, \ frac {7} {16} \ big) Überlegen. Das wahre $ p $ ist unbekannt, aber nur $ p_2 = p_3 $ ist bekannt. Zu diesem Zeitpunkt kann es ausgedrückt werden als $ p = [q, r, r, 1-2r-q] $. Finden Sie $ q, r $ mit der wahrscheinlichsten Schätzmethode.

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

Wird in einem blauen Histogramm angezeigt, und die theoretisch konvergierte $ \ chi ^ 2 $ -Verteilung mit einem Freiheitsgrad 1 wird in Orange angezeigt.

• n=4

• n=100

• n=10000

Klicken Sie hier für den Python-Quellcode, wenn $ 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)

Nächste Normalität der wahrscheinlichsten Schätzung (Teil 1)

Für die stochastische Variable $ X_i (i = 1, \ cdots, n) $, die durch den Parameter $ \ theta $ gekennzeichnet ist, ist $ \ theta_0 $ der wahre Parameter, $ \ hat \ theta $ ist die wahrscheinlichste Schätzung, Fisher-Information Betrag $ J_n (\ theta) $

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

Und. Wenn $ \ hat \ theta $ $ n \ bis \ infty $ ist, folgt es einer Normalverteilung des Mittelwerts $ \ theta_0 $ und der Varianz $ J_n (\ theta_0) ^ {-1} $.

Experiment

$ X_i (i = 1, \ cdots, n) $ Dichtefunktion

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

Und. Sei 10 das wahre $ \ theta $.

Die experimentell erhaltene $ \ hat \ theta $ -Verteilung ist blau und die theoretisch konvergierte Normalverteilung orange dargestellt.

• n=1

• n=10

• n=100

Klicken Sie hier, um den Python-Quellcode anzuzeigen, wenn $ n = 1 $.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
#Durchschnitt der Wahrscheinlichkeitsdichtefunktion 11/12,Dispersion 11/13-(11/12)^2
def f(x):
    return 11 * (x ** 10)
mu = 11 / 12
sigma2 = 11 / 13 - (11 / 12) ** 2
#Verteilungsfunktion
def F(x):
    return x ** 11
#Inverse Funktion der kumulativen Verteilungsfunktion
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)

#Theta-Histogramm
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]))
#Normalverteilung
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)

Nächste Normalität der wahrscheinlichsten Schätzung (Teil 2)

Sei $ g $ eine Funktion, die durch $ \ theta_0 $ unterschieden werden kann. Wenn $ g (\ hat \ theta) $ $ n \ bis \ infty $ ist, durchschnittlich $ g (\ theta_0) $, Dispersion $ g ^ \ prime (\ theta_0) ^ 2 J_n (\ theta_0) ^ {-1 } Folgen Sie der Normalverteilung von $. Die Definitionen von $ \ theta $, $ \ hat \ theta $ und $ \ theta_0 $ sind jedoch dieselben wie in Teil 1.

Experiment

Setzen Sie das Experiment von Teil 1 fort. Definieren Sie $ g $ mit der folgenden Formel.

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

Die experimentell erhaltene Verteilung von $ g (\ hat \ theta) $ ist blau und die theoretisch konvergierte Normalverteilung orange dargestellt.

• n=1

• n=10

• n=100

Klicken Sie hier für Python-Quellcode, wenn $ n = 1 $.

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

#Durchschnitt der Wahrscheinlichkeitsdichtefunktion 11/12,Dispersion 11/13-(11/12)^2
def f(x):
    return 11 * (x ** 10)
mu = 11 / 12
sigma2 = 11 / 13 - (11 / 12) ** 2
#Verteilungsfunktion
def F(x):
    return x ** 11
#Inverse Funktion der kumulativen Verteilungsfunktion
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)

Statistischer Test Betrag des Likelihood-Ratio-Tests

Angenommen, es gibt eine Wahrscheinlichkeitsverteilung $ f (x; \ theta) $, die durch den Parameter $ \ theta $ gekennzeichnet ist. Hypothesentestproblem

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

In wird das Wahrscheinlichkeitsverhältnis $ L $ durch die folgende Gleichung definiert.

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)}

Unter der Nullhypothese folgt $ 2 \ log L $ einer $ \ chi ^ 2 $ -Verteilung mit $ p $ Freiheit, wenn $ n \ bis \ infty $. $ P ＝ \ dim (\ Theta_1) - \ dim (\ Theta_0) $

Experiment

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

Als Hypothesentestproblem

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

Und.

Unter der Nullhypothese wird die Dichtefunktion von $ 2 \ log L $ im blauen Histogramm numerisch erhalten, und die Dichtefunktion der $ \ chi ^ 2 $ -Verteilung, die theoretisch konvergieren soll, ist die durchgezogene orange Linie.

• Wenn $ n = 1 $

• Wenn $ n = 10 $

Klicken Sie hier für den Quellcode, wenn $ n = 1 $.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import chi2
#Wahre Verteilung: Wahrscheinlichkeitsdichtefunktion Mittelwert 11/12,Dispersion 11/13-(11/12)^2
def f(x):
    return 11 * (x ** 10)

#Parametrisches Modell
def f(x, theta):
    return (theta + 1) * (x ** theta)

#Verteilungsfunktion
def F(x):
    return x ** 11

#Inverse Funktion der kumulativen Verteilungsfunktion
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  #Höchstwahrscheinlich Schätzung
    l.append(2 * np.log(np.prod(f(x, theta)) / np.prod(f(x, 10))))
l = np.array(l)
#l Histogramm
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]))
#Chi-Quadrat-Verteilung
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)