Quelles sont la moyenne et la variance d'une normale multivariée sans censure?

Nous pouvons d'abord réduire cela pour ne dépendre que de certains moments de distributions normales tronquées univariées / bivariées: notons bien sûr que $\DeclareMathOperator{\E}{\mathbb E} \DeclareMathOperator{\Var}{Var} \DeclareMathOperator{\Cov}{Cov} \newcommand{\N}{\mathcal N} \newcommand{\T}{\tilde} \newcommand{\v}{\mathcal V}$

\begin{matrix} E [Z_{+}] = {[\begin{matrix} E [(Z_{i})_{+}] \end{matrix}]}_{i} \\ Cov (Z_{+}) = {[\begin{matrix} Cov ((Z_{i})_{+}, (Z_{j})_{+}) \end{matrix}]}_{i j}, \end{matrix}

$\begin{gather} \E[Z_+] = \begin{bmatrix} \E[(Z_i)_+] \end{bmatrix}_i \\ \Cov(Z_+) = \begin{bmatrix} \Cov\left( (Z_i)_+, (Z_j)_+ \right) \end{bmatrix}_{ij} ,\end{gather}$ et parce que nous faisons des transformations coordonnées de certaines dimensions d'une distribution normale, nous ne besoin de s'inquiéter de la moyenne et de la variance d'une normale 1d censurée et de la covariance de deux normales 1d censurées.

Nous utiliserons quelques résultats de

S Rosenbaum (1961). Moments d'une distribution normale bivariée tronquée . JRSS B, vol 23 pp 405-408. ( jstor )

Rosenbaum considère et considère la troncature à l'événement .

[\begin{matrix} \tilde{X} \\ \tilde{Y} \end{matrix}] \sim N ([\begin{matrix} 0 \\ 0 \end{matrix}], [\begin{matrix} 1 & ρ \\ ρ & 1 \end{matrix}]),

$\begin{bmatrix}\T X \\ \T Y\end{bmatrix} \sim \N\left( \begin{bmatrix}0 \\ 0\end{bmatrix}, \begin{bmatrix}1 & \rho \\ \rho & 1\end{bmatrix} \right) ,$

V = {\tilde{X} \geq a_{X}, \tilde{Y} \geq a_{Y}}

$\v = \{ \T X \ge a_X, \T Y \ge a_Y \}$

Plus précisément, nous utiliserons les trois résultats suivants, ses (1), (3) et (5). Définissez d'abord les éléments suivants:

\begin{matrix} q_{x} = ϕ (a_{x}) q_{y} = ϕ (a_{y}) \\ Q_{x} = Φ (- a_{x}) Q_{y} = Φ (- a_{y}) \\ R_{x y} = Φ (\frac{ρ a_{x} - a_{y}}{\sqrt{1 - ρ^{2}}}) R_{y x} = Φ (\frac{ρ a_{y} - a_{x}}{\sqrt{1 - ρ^{2}}}) \\ r_{x y} = \frac{\sqrt{1 - ρ^{2}}}{\sqrt{2 π}} ϕ (\sqrt{\frac{h^{2} - 2 ρ h k + k^{2}}{1 - ρ^{2}}}) \end{matrix}

$\begin{gather} q_x = \phi( a_x) \qquad q_y = \phi( a_y) \\ Q_x = \Phi(-a_x) \qquad Q_y = \Phi(-a_y) \\ R_{xy} = \Phi\left( \frac{\rho a_x - a_y}{\sqrt{1 - \rho^2}} \right) \qquad R_{yx} = \Phi\left( \frac{\rho a_y - a_x}{\sqrt{1 - \rho^2}} \right) \\ r_{xy} = \frac{\sqrt{1-\rho^2}}{\sqrt{2 \pi}} \phi\left( \sqrt{\frac{h^2 - 2 \rho h k + k^2}{1 - \rho^2}} \right) \end{gather}$

Maintenant, Rosenbaum montre que:

\begin{aligned} (1) & Pr (V) E [\tilde{X} ∣ V] & = q_{x} R_{x y} + ρ q_{y} R_{y x} \\ (3) & Pr (V) E [{\tilde{X}}^{2} ∣ V] & = Pr (V) + a_{x} q_{x} R_{x y} + ρ^{2} a_{y} q_{y} R_{y x} + ρ r_{x y} \\ (5) & Pr (V) E [\tilde{X} \tilde{Y} ∣ V] & = ρ Pr (V) + ρ a_{x} q_{x} R_{x y} + ρ a_{y} q_{y} R_{y x} + r_{x y} . \end{aligned}

$\begin{align} \Pr(\v) \E[\T X \mid \v] &= q_x R_{xy} + \rho q_y R_{yx} \tag{1} \\ \Pr\left(\v \right) \E\left[\T X^2 \mid \v \right] &= \Pr\left(\v \right) + a_x q_x R_{xy} + \rho^2 a_y q_y R_{yx} + \rho r_{xy} \tag{3} \\ \Pr(\v) \E\left[ \T X \T Y \mid \v \right] &= \rho \Pr(\v) + \rho a_x q_x R_{xy} + \rho a_y q_y R_{yx} + r_{xy} \tag{5} .\end{align}$

Il sera utile de considérer également le cas particulier de (1) et (3) avec , c'est-à-dire une troncature 1d: $a_y = -\infty$

\begin{aligned} (*) & Pr (V) E [\tilde{X} ∣ V] & = q_{x} \\ (**) & Pr (V) E [{\tilde{X}}^{2} ∣ V] & = Pr (V) = Q_{x} . \end{aligned}

$\begin{align} \Pr(\v) \E[\T X \mid \v] &= q_x \tag{*} \\ \Pr\left(\v \right) \E\left[\T X^2 \mid \v \right] &= \Pr\left(\v \right) = Q_x \tag{**} .\end{align}$

Nous voulons maintenant considérer

\begin{aligned} [\begin{matrix} X \\ Y \end{matrix}] & = [\begin{matrix} μ_{x} \\ μ_{y} \end{matrix}] + [\begin{matrix} σ_{x} & 0 \\ 0 & σ_{y} \end{matrix}] [\begin{matrix} \tilde{X} \\ \tilde{Y} \end{matrix}] \\ \sim N ([\begin{matrix} μ_{X} \\ μ_{Y} \end{matrix}], [\begin{matrix} σ_{x}^{2} & ρ σ_{x} σ_{y} \\ ρ σ_{x} σ_{y} & σ_{y}^{2} \end{matrix}]) \\ = N (μ, Σ) . \end{aligned}

$\begin{align} \begin{bmatrix}X \\ Y\end{bmatrix} &= \begin{bmatrix}\mu_x\\\mu_y\end{bmatrix} + \begin{bmatrix}\sigma_x & 0 \\ 0 & \sigma_y\end{bmatrix}\begin{bmatrix}\T X \\ \T Y\end{bmatrix} \\&\sim \N\left( \begin{bmatrix}\mu_X\\\mu_Y\end{bmatrix}, \begin{bmatrix}\sigma_x^2 & \rho \sigma_x \sigma_y \\ \rho \sigma_x \sigma_y & \sigma_y^2 \end{bmatrix} \right) \\&= \N\left( \mu, \Sigma \right) .\end{align}$

Nous utiliserons qui sont les valeurs de et lorsque , .

a_{x} = \frac{- μ_{x}}{σ_{x}} a_{y} = \frac{- μ_{y}}{σ_{y}},

$a_x = \frac{-\mu_x}{\sigma_x} \qquad a_y = \frac{-\mu_y}{\sigma_y} ,$

\tilde{X}

$\T X$

\tilde{Y}

$\T Y$

X = 0

$X = 0$

Y = 0

$Y = 0$

Maintenant, en utilisant (*), nous obtenons et en utilisant à la fois (*) et (**) donne afin que

\begin{aligned} E [X_{+}] & = Pr (X_{+} > 0) E [X ∣ X > 0] + Pr (X_{+} = 0) 0 \\ = Pr (X > 0) (μ_{x} + σ_{x} E [\tilde{X} ∣ \tilde{X} \geq a_{x}]) \\ = Q_{x} μ_{x} + q_{x} σ_{x}, \end{aligned}

$\begin{align} \E[ X_+ ] &= \Pr(X_+ > 0) \E[X \mid X > 0] + \Pr(X_+=0) \, 0 \\&= \Pr(X > 0) \left( \mu_x + \sigma_x \E[\T X \mid \T X \ge a_x] \right) \\&= Q_x \mu_x + q_x \sigma_x ,\end{align}$

\begin{aligned} E [X_{+}^{2}] & = Pr (X_{+} > 0) E [X^{2} ∣ X > 0] + Pr (X_{+} = 0) 0 \\ = Pr (\tilde{X} \geq a_{x}) E [(μ_{x} + σ_{x} \tilde{X})^{2} ∣ \tilde{X} \geq a_{x}] \\ = Pr (\tilde{X} \geq a_{x}) E [μ_{x}^{2} + μ_{x} σ_{x} \tilde{X} + σ_{x}^{2} {\tilde{X}}^{2} ∣ \tilde{X} \geq a_{x}] \\ = Q_{x} μ_{x}^{2} + q_{x} μ_{x} σ_{x} + Q_{x} σ_{x}^{2} \end{aligned}

$\begin{align} \E[ X_+^2 ] &= \Pr(X_+ > 0) \E[X^2 \mid X > 0] + \Pr(X_+=0) 0 \\&= \Pr\left(\T X \ge a_x\right) \E\left[(\mu_x + \sigma_x \T X)^2 \mid \T X \ge a_x\right] \\&= \Pr\left(\T X \ge a_x\right) \E\left[\mu_x^2 + \mu_x \sigma_x \T X + \sigma_x^2 \T X^2 \mid \T X \ge a_x\right] \\&= Q_x \mu_x^2 + q_x \mu_x \sigma_x + Q_x \sigma_x^2 \end{align}$

\begin{aligned} Var [X_{+}] & = E [X_{+}^{2}] - E [X_{+}]^{2} \\ = Q_{x} μ_{x}^{2} + q_{x} μ_{x} σ_{x} + Q_{x} σ_{x}^{2} - Q_{x}^{2} μ_{x}^{2} - q_{x}^{2} σ_{x}^{2} - 2 q_{x} Q_{x} μ_{x} σ_{x} \\ = Q_{x} (1 - Q_{x}) μ_{x}^{2} + (1 - 2 Q_{x}) q_{x} μ_{x} σ_{x} + (Q_{x} - q_{x}^{2}) σ_{x}^{2} . \end{aligned}

$\begin{align} \Var[X_+] &= \E[X_+^2] - \E[X_+]^2 \\&= Q_x \mu_x^2 + q_x \mu_x \sigma_x + Q_x \sigma_x^2 - Q_x^2 \mu_x^2 - q_x^2 \sigma_x^2 - 2 q_x Q_x \mu_x \sigma_x \\&= Q_x (1 - Q_x) \mu_x^2 + (1 - 2 Q_x) q_x \mu_x \sigma_x + (Q_x - q_x^2) \sigma_x^2 .\end{align}$

Pour trouver , nous aurons besoin $\Cov(X_+, Y_+)$

\begin{aligned} E [X_{+} Y_{+}] & = Pr (V) E [X Y ∣ V] + P r (\neg V) 0 \\ = Pr (V) E [(μ_{x} + σ_{x} \tilde{X}) (μ_{y} + σ_{y} \tilde{Y}) ∣ V] \\ = μ_{x} μ_{y} Pr (V) + μ_{y} σ_{x} Pr (V) E [\tilde{X} ∣ V] + μ_{x} σ_{y} Pr (V) E [\tilde{Y} ∣ V] \\ + σ_{x} σ_{y} Pr (V) E [\tilde{X} \tilde{Y} ∣ V] \\ = μ_{x} μ_{y} Pr (V) + μ_{y} σ_{x} (q_{x} R_{x y} + ρ q_{y} R_{y x}) + μ_{x} σ_{y} (ρ q_{x} R_{x y} + q_{y} R_{y x}) \\ + σ_{x} σ_{y} (ρ Pr (V) - ρ μ_{x} q_{x} R_{x y} / σ_{x} - ρ μ_{y} q_{y} R_{y x} / σ_{y} + r_{x y}) \\ = (μ_{x} μ_{y} + σ_{x} σ_{y} ρ) Pr (V) + (μ_{y} σ_{x} + μ_{x} σ_{y} ρ - ρ μ_{x} σ_{y}) q_{x} R_{x y} \\ + (μ_{y} σ_{x} ρ + μ_{x} σ_{y} - ρ μ_{y} σ_{x}) q_{y} R_{y x} + σ_{x} σ_{y} r_{x y} \\ = (μ_{x} μ_{y} + Σ_{x y}) Pr (V) + μ_{y} σ_{x} q_{x} R_{x y} + μ_{x} σ_{y} q_{y} R_{y x} + σ_{x} σ_{y} r_{x y}, \end{aligned}

$\begin{align} \E[X_+ Y_+] &= \Pr(\v) \E[ X Y \mid \v] + Pr(\lnot\v) \, 0 \\&= \Pr(\v) \E\left[ (\mu_x + \sigma_x \T X) (\mu_y + \sigma_y \T Y) \mid \v \right] \\&= \mu_x \mu_y \Pr(\v) + \mu_y \sigma_x \Pr(\v) \E[ \T X \mid \v] + \mu_x \sigma_y \Pr(\v) \E[ \T Y \mid \v] \\&\qquad + \sigma_x \sigma_y \Pr(\v) \E\left[ \T X \T Y \mid \v \right] \\&= \mu_x \mu_y \Pr(\v) + \mu_y \sigma_x (q_x R_{xy} + \rho q_y R_{yx}) + \mu_x \sigma_y (\rho q_x R_{xy} + q_y R_{yx}) \\&\qquad + \sigma_x \sigma_y \left( \rho \Pr\left( \v \right) - \rho \mu_x q_x R_{xy} / \sigma_x - \rho \mu_y q_y R_{yx} / \sigma_y + r_{xy} \right) \\&= (\mu_x \mu_y + \sigma_x \sigma_y \rho) \Pr(\v) + (\mu_y \sigma_x + \mu_x \sigma_y \rho - \rho \mu_x \sigma_y) q_x R_{xy} \\&\qquad + (\mu_y \sigma_x \rho + \mu_x \sigma_y - \rho \mu_y \sigma_x) q_y R_{yx} + \sigma_x \sigma_y r_{xy} \\&= (\mu_x \mu_y + \Sigma_{xy}) \Pr(\v) + \mu_y \sigma_x q_x R_{xy} + \mu_x \sigma_y q_y R_{yx} + \sigma_x \sigma_y r_{xy} ,\end{align}$ puis en soustrayant nous obtenons

E [X_{+}] E [Y_{+}]

$\E[X_+] \E[Y_+]$

\begin{aligned} Cov (X_{+}, Y_{+}) & = (μ_{x} μ_{y} + Σ_{x y}) Pr (V) + μ_{y} σ_{x} q_{x} R_{x y} + μ_{x} σ_{y} q_{y} R_{y x} + σ_{x} σ_{y} r_{x y} \\ - (Q_{x} μ_{x} + q_{x} σ_{x}) (Q_{y} μ_{y} + q_{y} σ_{y}) . \end{aligned}

$\begin{align} \Cov(X_+, Y_+) &= (\mu_x \mu_y + \Sigma_{xy}) \Pr(\v) + \mu_y \sigma_x q_x R_{xy} + \mu_x \sigma_y q_y R_{yx} + \sigma_x \sigma_y r_{xy} \\&\qquad - (Q_x \mu_x + q_x \sigma_x) (Q_y \mu_y + q_y \sigma_y) .\end{align}$

Voici du code Python pour calculer les moments:

import numpy as np
from scipy import stats

def relu_mvn_mean_cov(mu, Sigma):
    mu = np.asarray(mu, dtype=float)
    Sigma = np.asarray(Sigma, dtype=float)
    d, = mu.shape
    assert Sigma.shape == (d, d)

    x = (slice(None), np.newaxis)
    y = (np.newaxis, slice(None))

    sigma2s = np.diagonal(Sigma)
    sigmas = np.sqrt(sigma2s)
    rhos = Sigma / sigmas[x] / sigmas[y]

    prob = np.empty((d, d))  # prob[i, j] = Pr(X_i > 0, X_j > 0)
    zero = np.zeros(d)
    for i in range(d):
        prob[i, i] = np.nan
        for j in range(i + 1, d):
            # Pr(X > 0) = Pr(-X < 0); X ~ N(mu, S) => -X ~ N(-mu, S)
            s = [i, j]
            prob[i, j] = prob[j, i] = stats.multivariate_normal.cdf(
                zero[s], mean=-mu[s], cov=Sigma[np.ix_(s, s)])

    mu_sigs = mu / sigmas

    Q = stats.norm.cdf(mu_sigs)
    q = stats.norm.pdf(mu_sigs)
    mean = Q * mu + q * sigmas

    # rho_cs is sqrt(1 - rhos**2); but don't calculate diagonal, because
    # it'll just be zero and we're dividing by it (but not using result)
    # use inf instead of nan; stats.norm.cdf doesn't like nan inputs
    rho_cs = 1 - rhos**2
    np.fill_diagonal(rho_cs, np.inf)
    np.sqrt(rho_cs, out=rho_cs)

    R = stats.norm.cdf((mu_sigs[y] - rhos * mu_sigs[x]) / rho_cs)

    mu_sigs_sq = mu_sigs ** 2
    r_num = mu_sigs_sq[x] + mu_sigs_sq[y] - 2 * rhos * mu_sigs[x] * mu_sigs[y]
    np.fill_diagonal(r_num, 1)  # don't want slightly negative numerator here
    r = rho_cs / np.sqrt(2 * np.pi) * stats.norm.pdf(np.sqrt(r_num) / rho_cs)

    bit = mu[y] * sigmas[x] * q[x] * R
    cov = (
        (mu[x] * mu[y] + Sigma) * prob
        + bit + bit.T
        + sigmas[x] * sigmas[y] * r
        - mean[x] * mean[y])

    cov[range(d), range(d)] = (
        Q * (1 - Q) * mu**2 + (1 - 2 * Q) * q * mu * sigmas
        + (Q - q**2) * sigma2s)

    return mean, cov

et un test de Monte Carlo que cela fonctionne:

np.random.seed(12)
d = 4
mu = np.random.randn(d)
L = np.random.randn(d, d)
Sigma = L.T.dot(L)
dist = stats.multivariate_normal(mu, Sigma)

mn, cov = relu_mvn_mean_cov(mu, Sigma)

samps = dist.rvs(10**7)
mn_est = samps.mean(axis=0)
cov_est = np.cov(samps, rowvar=False)
print(np.max(np.abs(mn - mn_est)), np.max(np.abs(cov - cov_est)))

ce qui donne 0.000572145310512 0.00298692620286, indiquant que l'espérance et la covariance revendiquées correspondent aux estimations de Monte Carlo (sur la base de échantillons). $10,000,000$

Dougal
la source

pouvez-vous résumer quelles sont ces valeurs finales? S'agit-il d'estimations des paramètres mu et L que vous avez générés? Peut-être imprimer ces valeurs cibles?

AdamO

Non, les valeurs de retour sont et ; ce que j'ai imprimé était la distance entre les estimateurs de Monte Carlo de ces quantités et la valeur calculée. Vous pouvez peut-être inverser ces expressions pour obtenir un estimateur de correspondance de moment pour et - Rosenbaum le fait en fait dans sa section 3 dans le cas tronqué - mais ce n'est pas ce que je voulais ici.

\E (Z_{+})

$\E(Z_+)$

\Cov (Z_{+})

$\Cov(Z_+)$

L_{\infty}

$L_\infty$

μ

$\mu$

Σ

$\Sigma$

Dougal

Quelles sont la moyenne et la variance d'une normale multivariée sans censure?

Réponses: