Dérivation de la distribution t multivariée et implémentation de la génération de nombres aléatoires par python

introduction

Puisque j'étudiais la distribution multivariée t, c'est un résumé.

(Au moment où j'ai écrit cet article, j'ai soufflé les données et je les ai réécrites ... je veux mourir ...)

Dérivation de la distribution t multivariée

Définition

Suivez la notation de Distribution t multivariée sur wikipedia car Gugu est le premier à être trouvé. La fonction de densité de probabilité pour la distribution t multivariée est définie ci-dessous.

théorie

point

En premier lieu, qu'est-ce que la distribution t? On peut dire que la distribution t est un paramètre qui augmente ou diminue la dispersion de la distribution normale, suppose une distribution gamma et finalement élimine l'intégration. Ici, nous pouvons supposer la distribution gamma inverse et la distribution du chi carré au lieu de la distribution gamma, mais il semble qu'elles auront toutes la même valeur avec un léger changement. (Je n'ai pas essayé)

Cela vous permet de voir le comportement moyen lorsque vous n'êtes pas sûr de la valeur de la variance.

Dérivation

Considérons deux variables stochastiques $ V \ sim Gam (\ alpha, \ beta) $, $ X \ sim N (\ textbf {0}, {\ bf \ Sigma}) $. $ X $ est une distribution normale de dimension $ d $. Ces deux variables stochastiques sont converties en $ Y et U $ par les fonctions $ f et g . $ Y = f(V,X) = \mu + \frac{1}{\sqrt{V}}X, \tag{2} $ $ U = g(V,X) = g(V)= V. \tag{3} $$

Ici, $ g $ est une fonction égale de $ V $, $ \ mu \ in \ mathbb {R} ^ {d} $ est un vecteur à d dimensions représentant la moyenne, et $ V $ est indépendant de $ X $. Je vais assumer.

Cette $ f (\ cdot) $ est une fonction de ligne droite (plan) lorsque $ V $ est fixe, et est une fonction de correspondance un-à-un. Par conséquent, la fonction inverse $ f ^ {-1} (\ cdot) $ existe, $ f^{-1}(Y) = X = \sqrt{U}(Y - \mu), \tag{4} $ Il devient.

Trouvez la distribution de $ Y $ à partir de $ V et X $ dont la distribution est déjà connue. Puisqu'il existe deux variables stochastiques, $ U $ et $ Y $, la distribution de $ Y $ est obtenue en intégrant et en éliminant $ U $ en considérant la distribution simultanée de $ U et Y $.

Pr(u,y) = Pr(U < u, Y < y) = Pr(g(V) < u,f(X) < y) \\
= Pr(V < g^{-1}(u), X < f^{-1}(y)) \\
= Pr(V < u, X<\sqrt{u}(y - \mu)) = Pr(V < u) Pr(X < \sqrt{u}(y - \mu)). \tag{5}

Ici, $ V $ et $ X $ sont supposés indépendants, $ Pr(v,x) = Pr(v)Pr(x), \tag{6} $ J'ai utilisé.

Soit $ p_ {V} et p_ {X} $ les fonctions de densité de probabilité de $ V $ et $ X , respectivement. $ Pr(u,y) = Pr(V < u) Pr(X < \sqrt{u}(y - \mu)) = \int_{-\infty}^{u} \int_{-\infty}^{y} p_{V}(u)p_{X} (\sqrt{u}(y - \mu)) |\det(J)| dudy. \tag{7} $$

ici,|\det(J)|Le jacobien est-il généré par transformation variable et la matrice de JacobiJ （(d+1)*(d+1)Matrice)

J = \begin{pmatrix} \frac{\partial x}{\partial y^{T}} & \frac{\partial x}{\partial u} \\ \frac{\partial v}{\partial y^{T}} & \frac{\partial v}{\partial u} \end{pmatrix} =
 \begin{pmatrix} \sqrt{u}I_{d} & \frac{1}{2\sqrt{u}(y-\mu)} \\ \textbf{0} & 1 \end{pmatrix}. \tag{8}

Puisque la matrice de Jacobi $ J $ est une matrice triangulaire supérieure, l'expression matricielle est le produit d'éléments diagonaux. $ \det(J) = u^{\frac{d}{2}}. \tag{9} $ 　　 Puisque le contenu de l'intégration de $ (7) $ est la fonction de densité de probabilité $ p_ {U, Y} (u, y) $ de $ Pr (u, y) $, ceci est calculé. À partir de $ V \ sim Gam (\ alpha, \ beta) $, $ X \ sim N (\ textbf {0}, {\ bf \ Sigma}) $ $ p_{V}(v) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} v^{\alpha-1} \exp{( - \beta v)}, \tag{10} $ $ p_{X}(\textbf{x}) = \frac{1}{(2\pi)^{\frac{d}{2}}} \frac{1}{\sqrt{\det(\Sigma)}} \exp{\Bigl( - \frac{1}{2} \textbf{x}^{T} \Sigma^{-1} \textbf{x} \Bigr)}, \tag{11} $ Donc, $ p_{V}(u) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} u^{\alpha-1} \exp{( - \beta u)} = Gam(\alpha, \beta), \tag{12} $

p_{X}(\sqrt{u}(\textbf{y} - \mu)) = \frac{1}{(2\pi)^{\frac{d}{2}}} \frac{1}{\sqrt{\det(\Sigma)}} \exp{\Bigl( - \frac{1}{2} \sqrt{u}(\textbf{y} - \mu)^{T} \Sigma^{-1} \sqrt{u}(\textbf{y} - \mu) \Bigr)} \\
= \frac{1}{(2\pi)^{\frac{d}{2}}} u^{-\frac{d}{2}} u^{\frac{d}{2}} \frac{1}{\sqrt{\det(\Sigma)}} \exp{\Bigl( - \frac{1}{2} (\textbf{y} - \mu)^{T} (u \Sigma^{-1}) (\textbf{y} - \mu) \Bigr)} \\
=  u^{-\frac{d}{2}} \frac{1}{(2\pi)^{\frac{d}{2}}}  \frac{1}{\sqrt{\det(\frac{\Sigma}{u})}} \exp{\Bigl( - \frac{1}{2} (\textbf{y} - \mu)^{T} (u \Sigma^{-1}) (\textbf{y} - \mu) \Bigr)} = u^{-\frac{d}{2}} N(\mu, \frac{\Sigma}{u}). \tag{13}

La fonction de densité de probabilité $ p_ {U, Y} (u, y) $ de $ Pr (u, y) $ est $ p_{U,Y}(u,y) = Gam(\alpha, \beta) u^{-\frac{d}{2}} N(\mu, \frac{\Sigma}{u}) |\det(J)| = Gam(\alpha, \beta) N(\mu, \frac{\Sigma}{u}) u^{-\frac{d}{2}}u^{\frac{d}{2}} = Gam(\alpha, \beta) N(\mu, \frac{\Sigma}{u}). \tag{14} $ (Ici, la plateforme de la distribution gamma est $ u \ in [0, \ infty), d> 0 $, donc la valeur absolue peut être supprimée. )

En intégrant $ p_ {U, Y} (u, y) $ par rapport à $ u $, on obtient la fonction de densité de probabilité $ p_ {Y} $ de $ Y . $ p_{Y} = \int_{0}^{\infty} p_{U,Y}(u,y)du = \int_{0}^{\infty} Gam(\alpha, \beta) N(\mu, \frac{\Sigma}{u})du \tag{15} $$

Ici, par souci de simplification de la notation, la distance du Maharanobis normalement distribué est fixée à $ \ Delta ^ {2} = (\ textbf {y} - \ mu) ^ {T} \ Sigma ^ {-1} (). Il s'écrit textbf {y} - \ mu) $.

Extraire uniquement la partie qui dépend de $ u $ $ p_{Y}(\textbf{y}) \propto \int_{0}^{\infty} u^{\alpha-1} \exp{( - \beta u)} u^{\frac{d}{2}} \exp{\Bigl( - \frac{u}{2} \Delta^{2} \Bigr)} du \\ = \int_{0}^{\infty} u^{\frac{d}{2}+\alpha-1} \exp{\Bigl( - u( \beta + \frac{\Delta^{2}}{2} ) \Bigr)} du. \tag{17} $

Ici, si $ z = u (\ beta + \ frac {\ Delta ^ {2}} {2}) $, alors $ u = z (\ beta + \ frac {\ Delta ^ {2}} {2 }) ^ {-1} $, et $ du = (\ beta + \ frac {\ Delta ^ {2}} {2}) ^ {-1} dz $. Lorsque $ u = 0 $, $ z = 0 $, $ u \ to \ infty $ et $ z \ to \ infty $.

De la définition de la fonction gamma $ p_{Y}(\textbf{y}) \propto \Bigl( \beta + \frac{\Delta^{2}}{2} \Bigr)^{-\frac{d}{2} -\alpha} \Gamma \Bigl( \frac{d}{2} + \alpha \Bigr) \tag{18} $

Compte tenu du terme indépendant de $ u $ de (17) $

p_{Y}(\textbf{y}) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} \frac{1}{(2\pi)^{\frac{d}{2}}}  \frac{1}{\sqrt{\det(\Sigma)}} \Bigl( \beta + \frac{\Delta^{2}}{2} \Bigr)^{-\frac{d}{2} -\alpha}  \\
= \frac{\Gamma \Bigl( \frac{d}{2} + \alpha \Bigr)}{\Gamma(\alpha)}\frac{\beta^{\alpha}}{(2\pi)^{\frac{d}{2}}} \frac{1}{\sqrt{\det(\Sigma)}} \Bigl( \beta + \frac{\Delta^{2}}{2} \Bigr)^{-\frac{d}{2} -\alpha} \Bigl( \frac{1}{\beta} \Bigr)^{\frac{d}{2} + \alpha} \beta^{-\frac{d}{2}-\alpha} \\
= \frac{\Gamma \Bigl( \frac{d}{2} + \alpha \Bigr)}{\Gamma(\alpha)}\frac{1}{(2 \beta \pi)^{\frac{d}{2}}} \frac{1}{\sqrt{\det(\Sigma)}} \Bigl( 1 + \frac{\Delta^{2}}{2 \beta} \Bigr)^{-\frac{d}{2} -\alpha}. \tag{19}

Ici, si $ \ alpha = \ beta = \ frac {\ nu} {2} $ et $ p = d $ sont définis, $ p_{Y}(\textbf{y}) = \frac{\Gamma \Bigl( \frac{\nu+p}{2} \Bigr)}{\Gamma \Bigl( \frac{\nu}{2} \Bigr)}\frac{1}{(\nu \pi)^{\frac{p}{2}}} \frac{1}{\sqrt{\det(\Sigma)}} \Bigl( 1 + \frac{\Delta^{2}}{\nu} \Bigr)^{-\frac{(\nu+p)}{2}}, \tag{20} $ Est en accord avec la définition de la distribution t de $ (1) $. $ \ Boîte $

Entraine toi

Implémentation de la génération de nombres aléatoires par Python

Puisque l'expression $ (2) $ est la méthode de génération de nombres aléatoires telle qu'elle est, nous allons l'implémenter en python.

# -*- coding: utf-8 -*-
#!/usr/bin/python

import numpy as np
import matplotlib.pyplot as plt

#Réglez les paramètres.
N = 5000
d = 10
df = 4

mean = np.zeros(d)
cov = np.eye(d)

#Génération aléatoire à partir de la distribution gamma
V = np.random.gamma(shape=df/2., scale=df/2.,size=N)

#Génération aléatoire à partir d'une distribution normale multivariée
X = np.random.multivariate_normal(mean=mean, cov=cov, size=N)

#Convertissez X.
V = V.reshape((N,1))
denom = np.sqrt(np.tile(V,d))
Y = mean + X / denom


#Générez des nombres aléatoires à partir d'une distribution normale multivariée pour comparaison.
#(Question:Est-il acceptable que la comparaison ressemble à ceci? )
cov_ = ((df - 2) / float(df)) * np.corrcoef(Y.T)
X_ = np.random.multivariate_normal(mean=mean, cov=cov_, size=N)

#Terrain.
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hist(X_[:,0], color="blue", normed=True, label="Normal dist", bins=30)
ax.hist(Y[:,0], color="red", normed=True, label="t dist", bins=30, alpha=0.5)
ax.set_title("Comparison of normal and t distribution")
ax.legend()
fig.show()

fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(X_[:,0], X_[:,1],color="blue",label="Normal dist")
ax.scatter(Y[:,0], Y[:,1],color="red",label="t dist",alpha=0.5)
ax.set_title("scatter plot of t and normal dist")
ax.legend()
fig.show()

La dispersion est grande et cela ressemble à une distribution t? (Je voulais comparer autre chose qu'une dimension à la fois, mais je l'ai réécrit et épuisé, donc c'est à peu près tout ...)

finalement

Atmosphere Mathematics man, alors faites le moi savoir si vous faites une erreur. Si j'ai le temps, j'aimerais suivre la connexion des différentes distributions.