Was sind der Mittelwert und die Varianz einer 0-zensierten multivariaten Normalen?

9

ZN(μ,Σ)RdZ+=max(0,Z)

Dies tritt z. B. auf, wenn wir die ReLU-Aktivierungsfunktion in einem tiefen Netzwerk verwenden und über das CLT annehmen, dass die Eingaben in eine bestimmte Schicht ungefähr normal sind, dann ist dies die Verteilung der Ausgaben.

(Ich bin sicher, dass viele Leute dies bereits berechnet haben, aber ich konnte das aufgelistete Ergebnis nirgends auf eine einigermaßen lesbare Weise finden.)

Dougal
quelle
Es würde Ihre Antwort - vielleicht sehr - vereinfachen, zu beobachten, dass Sie sie erhalten können, indem Sie die Ergebnisse zweier getrennter Fragen kombinieren: (1) Was sind die Momente einer abgeschnittenen Normalverteilung und (2) Was sind die Momente einer Mischung? ? Letzteres ist unkompliziert und alles, was Sie tun müssen, ist, Ergebnisse für Ersteres zu zitieren.
whuber
@whuber Hmm. Obwohl ich es nicht explizit gesagt habe, ist es im Wesentlichen das, was ich in meiner Antwort mache, außer dass ich keine Ergebnisse für eine abgeschnittene bivariate Verteilung mit einem allgemeinen Mittelwert und einer allgemeinen Varianz gefunden habe und deshalb trotzdem etwas skalieren und verschieben musste. Gibt es eine Möglichkeit, z. B. die Kovarianz abzuleiten, ohne die Menge an Algebra zu tun, die ich tun musste? Ich behaupte mit Sicherheit nicht, dass irgendetwas in dieser Antwort neu ist, nur dass die Algebra langweilig und fehleranfällig war und vielleicht jemand anderes die Lösung nützlich finden wird.
Dougal
Richtig: Ich bin sicher, dass Ihre Algebra gleichbedeutend mit dem ist, was ich beschrieben habe. Es sieht also so aus, als ob wir die Wertschätzung teilen, dass eine Vereinfachung der Algebra möglich ist. Eine einfache Möglichkeit, die Algebra zu reduzieren, besteht darin, die diagonalen Elemente von auf Eins zu standardisieren , da lediglich für jede Variable eine Maßeinheit festgelegt wird. An diesem Punkt können Sie Rosenbaums Ergebnisse direkt in die (einfachen, offensichtlichen) Ausdrücke für Momente von Gemischen einfügen. Ob dies überhaupt eine algebraische Vereinfachung wert ist, mag Geschmackssache sein: Ohne Vereinfachung führt dies zu einem einfachen, modularen Computerprogramm. Σ
whuber
1
Ich nehme an, man könnte ein Programm schreiben, das Momente direkt mit Rosenbaums Ergebnissen berechnet und angemessen mischt und sie dann verschiebt und zurück in den ursprünglichen Raum skaliert. Das wäre wahrscheinlich schneller gewesen als ich es getan habe.
Dougal

Antworten:

7

Wir können dies zunächst so reduzieren, dass es nur von bestimmten Momenten univariater / bivariater abgeschnittener Normalverteilungen abhängt: Beachten Sie natürlich das

E[Z+]=[E[(Zi)+]]iCov(Z+)=[Cov((Zi)+,(Zj)+)]ij,
und weil wir koordinatweise Transformationen bestimmter Dimensionen einer Normalverteilung durchführen, nur wir Sie müssen sich um den Mittelwert und die Varianz einer 1d-zensierten Normalen und die Kovarianz zweier 1d-zensierter Normalen sorgen.

Wir werden einige Ergebnisse von verwenden

S. Rosenbaum (1961). Momente einer abgeschnittenen bivariaten Normalverteilung . JRSS B, Bd. 23, S. 405-408. ( jstor )

Rosenbaum betrachtet und berücksichtigt das Abschneiden auf das Ereignis . V ={ ˜ X a X , ˜ Y a Y }

[X~Y~]N([00],[1ρρ1]),
V={X~aX,Y~aY}

Insbesondere werden wir die folgenden drei Ergebnisse verwenden, seine (1), (3) und (5). Definieren Sie zunächst Folgendes:

qx=ϕ(ax)qy=ϕ(ay)Qx=Φ(ax)Qy=Φ(ay)Rxy=Φ(ρaxay1ρ2)Ryx=Φ(ρayax1ρ2)rxy=1ρ22πϕ(h22ρhk+k21ρ2)

Nun zeigt Rosenbaum, dass:

(1)Pr(V)E[X~V]=qxRxy+ρqyRyx(3)Pr(V)E[X~2V]=Pr(V)+axqxRxy+ρ2ayqyRyx+ρrxy(5)Pr(V)E[X~Y~V]=ρPr(V)+ρaxqxRxy+ρayqyRyx+rxy.

Es ist nützlich, auch den Sonderfall von (1) und (3) mit , dh eine 1d-Kürzung: Pr ( V ) E [ ˜ XV ]ay=

(*)Pr(V)E[X~V]=qx(**)Pr(V)E[X~2V]=Pr(V)=Qx.

Wir wollen nun

[XY]=[μxμy]+[σx00σy][X~Y~]N([μXμY],[σx2ρσxσyρσxσyσy2])=N(μ,Σ).

Wir werden sind die Werte von und wenn , .

ax=μxσxay=μyσy,
X~Y~X=0Y=0

Mit (*) erhalten wir nun und die Verwendung von (*) und (**) ergibt

E[X+]=Pr(X+>0)E[XX>0]+Pr(X+=0)0=Pr(X>0)(μx+σxE[X~X~ax])=Qxμx+qxσx,
E[X+2]=Pr(X+>0)E[X2X>0]+Pr(X+=0)0=Pr(X~ax)E[(μx+σxX~)2X~ax]=Pr(X~ax)E[μx2+μxσxX~+σx2X~2X~ax]=Qxμx2+qxμxσx+Qxσx2
so dass
Var[X+]=E[X+2]E[X+]2=Qxμx2+qxμxσx+Qxσx2Qx2μx2qx2σx22qxQxμxσx=Qx(1Qx)μx2+(12Qx)qxμxσx+(Qxqx2)σx2.

Um , benötigen wir Cov(X+,Y+)

E[X+Y+]=Pr(V)E[XYV]+Pr(¬V)0=Pr(V)E[(μx+σxX~)(μy+σyY~)V]=μxμyPr(V)+μyσxPr(V)E[X~V]+μxσyPr(V)E[Y~V]+σxσyPr(V)E[X~Y~V]=μxμyPr(V)+μyσx(qxRxy+ρqyRyx)+μxσy(ρqxRxy+qyRyx)+σxσy(ρPr(V)ρμxqxRxy/σxρμyqyRyx/σy+rxy)=(μxμy+σxσyρ)Pr(V)+(μyσx+μxσyρρμxσy)qxRxy+(μyσxρ+μxσyρμyσx)qyRyx+σxσyrxy=(μxμy+Σxy)Pr(V)+μyσxqxRxy+μxσyqyRyx+σxσyrxy,
und dann subtrahiert man so erhält man E[X+]E[Y+]
Cov(X+,Y+)=(μxμy+Σxy)Pr(V)+μyσxqxRxy+μxσyqyRyx+σxσyrxy(Qxμx+qxσx)(Qyμy+qyσy).

Hier ist ein Python-Code, um die Momente zu berechnen:

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

und ein Monte-Carlo-Test, dass es funktioniert:

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

Dies gibt an 0.000572145310512 0.00298692620286, dass die behauptete Erwartung und Kovarianz mit den Monte-Carlo-Schätzungen übereinstimmen (basierend auf Stichproben).10,000,000

Dougal
quelle
Können Sie zusammenfassen, was diese Endwerte sind? Sind sie Schätzungen der von Ihnen generierten Parameter mu und L? Vielleicht diese Zielwerte drucken?
AdamO
Nein, die Rückgabewerte sind und ; Was ich war der gedruckte Abstand zwischen Monte Carlo Schätzer dieser Mengen und dem berechneten Wert. Sie könnten diese Ausdrücke vielleicht invertieren, um einen Schätzer für die Momentanpassung für und - Rosenbaum tut dies tatsächlich in seinem Abschnitt 3 im abgeschnittenen Fall -, aber das wollte ich hier nicht. \ Cov ( Z + ) L μ Σ\E(Z+)\Cov(Z+)LμΣ
Dougal