Jaka jest średnia i wariancja normalnej wielowymiarowej 0-cenzurowanej?

9

Niech będzie w . Jaka jest średnia i macierz kowariancji (z maksimum obliczonym elementarnie)?ZN(μ,Σ)RdZ+=max(0,Z)

Dzieje się tak np. Dlatego, że jeśli użyjemy funkcji aktywacji ReLU w głębokiej sieci i założymy przez CLT, że wejścia do danej warstwy są w przybliżeniu normalne, to jest to rozkład wyjść.

(Jestem pewien, że wiele osób już to obliczało, ale nie mogłem znaleźć nigdzie wymienionego wyniku w sposób czytelny).

Dougal
źródło
Uprościłoby to twoją odpowiedź - być może znacznie - obserwując, że możesz ją uzyskać, łącząc wyniki dwóch oddzielnych pytań: (1) jakie są momenty okrojonego rozkładu normalnego i (2) jakie są momenty mieszanki ? To drugie jest proste i wszystko, co musisz zrobić, to przytoczyć wyniki pierwszego.
whuber
@whuber Hmm. Chociaż nie powiedziałem tego wprost, to zasadniczo robię w mojej odpowiedzi, z tym wyjątkiem, że nie znalazłem wyników dla skróconego rozkładu dwuwymiarowego z ogólną średnią i wariancją, więc i tak musiałem trochę skalować i przesuwać. Czy jest jakiś sposób na uzyskanie np. Kowariancji bez wykonywania algebry, którą musiałem zrobić? Z pewnością nie twierdzę, że cokolwiek w tej odpowiedzi jest nowe, tylko że algebra była żmudna i podatna na błędy i być może ktoś inny uzna to rozwiązanie za przydatne.
Dougal
Racja: Jestem pewien, że twoja algebra jest równoznaczna z tym, co opisałem, więc wygląda na to, że podzielamy uznanie dla uproszczenia algebry. Jednym łatwym sposobem na zmniejszenie algebry jest standaryzacja elementów ukośnych do jedności, ponieważ wystarczy ustalić jednostkę miary dla każdej zmiennej. W tym momencie możesz bezpośrednio podłączyć wyniki Rosenbauma do (prostych, oczywistych) wyrażeń dla momentów mieszanin. To, czy warto w ogóle algebraiczne uproszczenie, mogło być kwestią gustu: bez uproszczenia prowadzi do prostego, modułowego programu komputerowego. Σ
whuber
1
Przypuszczam, że można napisać program, który oblicza momenty bezpośrednio z wynikami Rosenbauma i odpowiednio je miksuje, a następnie przesuwa je i skaluje z powrotem do pierwotnej przestrzeni. To pewnie byłoby szybsze niż ja to zrobiłem.
Dougal

Odpowiedzi:

7

Możemy to najpierw zmniejszyć, aby zależeć tylko od pewnych momentów jednostronnego / dwuwymiarowego skróconego rozkładu normalnego: pamiętaj oczywiście, że

E[Z+]=[E[(Zi)+]]iCov(Z+)=[Cov((Zi)+,(Zj)+)]ij,
a ponieważ dokonujemy transformacji współrzędnych niektórych wymiarów rozkładu normalnego, tylko trzeba się martwić o średnią i wariancję normalnej ocenzurowanej 1d oraz kowariancji dwóch normalnych ocenzurowanych 1d.

Wykorzystamy niektóre wyniki z

S Rosenbaum (1961). Momenty skróconego dwuwymiarowego rozkładu normalnego . JRSS B, tom 23 str. 405–408. ( jstor )

Rosenbaum uważa, że i rozważa obcięcie zdarzenia .

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

W szczególności użyjemy trzech następujących wyników, jego (1), (3) i (5). Najpierw zdefiniuj następujące elementy:

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

Teraz Rosenbaum pokazuje, że:

(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.

Przydatne będzie również rozważenie specjalnego przypadku (1) i (3) z , tj. Obcinaniem 1d: ay=

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

Teraz chcemy rozważyć

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

Użyjemy które są wartościami i gdy , .

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

Teraz, używając (*), otrzymujemy i użycie zarówno (*), jak i (**) daje dzięki czemu

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

Aby znaleźć , będziemy potrzebować 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,
a następnie odejmując otrzymujemy 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).

Oto kod Pythona do obliczenia chwil:

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

oraz test Monte Carlo, że działa:

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

co daje 0.000572145310512 0.00298692620286, wskazując, że deklarowane oczekiwania i kowariancja odpowiadają szacunkom Monte Carlo (na podstawie próbek).10,000,000

Dougal
źródło
czy możesz podsumować, jakie są te końcowe wartości? Czy są to oszacowania parametrów mu i L, które wygenerowałeś? Może wydrukujesz te wartości docelowe?
AdamO,
Nie, zwracanymi wartościami są i ; wydrukowałem odległość między estymatorami Monte Carlo tych ilości a wartością obliczoną. Możesz odwrócić te wyrażenia, aby uzyskać estymator dopasowujący moment dla i - Rosenbaum faktycznie robi to w swojej sekcji 3 w skróconej sprawie - ale nie tego tu chciałem. \E(Z+)\Cov(Z+)LμΣ
Dougal