Apa arti dan varian dari normal multivariat 0-disensor?

9

Biarkan berada di . Apa yang dimaksud dengan matriks kovarians dan rata-rata dari (dengan max dihitung dengan elementwise)?ZN(μ,Σ)RdZ+=max(0,Z)

Ini muncul misalnya karena, jika kita menggunakan fungsi aktivasi ReLU di dalam jaringan yang dalam, dan mengasumsikan melalui CLT bahwa input ke lapisan yang diberikan kira-kira normal, maka ini adalah distribusi dari output.

(Saya yakin banyak orang telah menghitung ini sebelumnya, tetapi saya tidak dapat menemukan hasil yang tercantum di mana pun dengan cara yang cukup mudah dibaca.)

Dougal
sumber
Ini akan menyederhanakan jawaban Anda - mungkin sangat - untuk mengamati bahwa Anda dapat memperolehnya dengan menggabungkan hasil dari dua pertanyaan terpisah: (1) apa saja momen dari distribusi Normal terpotong dan (2) apa saat-saat campuran ? Yang terakhir ini mudah dan yang perlu Anda lakukan adalah mengutip hasil untuk yang pertama.
whuber
@whuber Hmm. Meskipun saya tidak mengatakannya secara eksplisit, itu pada dasarnya apa yang saya lakukan dalam jawaban saya, kecuali bahwa saya tidak menemukan hasil untuk distribusi bivariat terpotong dengan rata-rata dan varian umum sehingga harus melakukan penskalaan dan penggeseran pula. Apakah ada cara untuk menurunkan misalnya kovarians tanpa melakukan jumlah aljabar yang harus saya lakukan? Saya tentu tidak mengklaim bahwa apa pun dalam jawaban ini adalah baru, hanya saja aljabarnya membosankan dan rentan kesalahan, dan mungkin orang lain akan menemukan solusi yang bermanfaat.
Dougal
Kanan: Saya yakin aljabar Anda sama dengan apa yang saya jelaskan, jadi sepertinya kami berbagi penghargaan untuk menyederhanakan aljabar itu mungkin. Salah satu cara mudah untuk mengurangi aljabar adalah dengan membakukan elemen diagonal menjadi satu, karena semua yang dilakukan adalah membuat unit pengukuran untuk setiap variabel. Pada titik itu Anda bisa langsung menyambungkan hasil Rosenbaum ke dalam ekspresi (sederhana, jelas) untuk momen campuran. Apakah itu berarti penyederhanaan aljabar yang layak mungkin hanya masalah selera: tanpa penyederhanaan, ini mengarah pada program komputer modular yang sederhana. Σ
whuber
1
Saya kira orang bisa menulis sebuah program yang menghitung momen secara langsung dengan hasil Rosenbaum dan pencampuran dengan tepat, dan kemudian menggeser dan menskalakannya kembali ke ruang asli. Itu mungkin akan lebih cepat daripada cara saya melakukannya.
Dougal

Jawaban:

7

Pertama-tama kita dapat mengurangi ini hanya bergantung pada momen-momen tertentu dari distribusi normal univariat / bivariat terpotong: perhatikan tentu saja bahwa

E[Z+]=[E[(Zi)+]]iCov(Z+)=[Cov((Zi)+,(Zj)+)]ij,
dan karena kami membuat transformasi bijaksana dimensi tertentu dari distribusi normal, kami hanya perlu khawatir tentang mean dan varians dari sensor normal 1d disensor dan kovarians dari dua sensor normal 1d disensor.

Kami akan menggunakan beberapa hasil dari

S Rosenbaum (1961). Momen dari Distribusi Normal Bivariat Terpotong . JRSS B, vol 23 hlm 405-408. ( jstor )

Rosenbaum menganggap dan pertimbangkan pemotongan ke acara .

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

Secara khusus, kami akan menggunakan tiga hasil berikut ini, (1), (3), dan (5). Pertama, tentukan yang berikut:

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

Sekarang, Rosenbaum menunjukkan bahwa:

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

Akan berguna juga untuk mempertimbangkan kasus khusus (1) dan (3) dengan , yaitu pemotongan 1d: ay=

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

Kami sekarang ingin mempertimbangkan

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

Kita akan menggunakan yang merupakan nilai dan ketika , .

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

Sekarang, menggunakan (*), kita memperoleh dan menggunakan keduanya (*) dan (**) menghasilkan sehingga

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.

Untuk menemukan , kita perlu 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,
dan kemudian mengurangi kita mendapatkan 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).

Berikut ini beberapa kode Python untuk menghitung momen:

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

dan tes Monte Carlo yang berfungsi:

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

yang memberi 0.000572145310512 0.00298692620286, menunjukkan bahwa ekspektasi dan kovarian yang diklaim sesuai dengan perkiraan Monte Carlo (berdasarkan sampel).10,000,000

Dougal
sumber
dapatkah Anda meringkas apa nilai akhir itu? Apakah mereka memperkirakan parameter mu dan L yang Anda hasilkan? Mungkin mencetak nilai target itu?
AdamO
Tidak, nilai yang dikembalikan adalah dan ; apa yang saya cetak adalah jarak tak antara estimasi Monte Carlo dari jumlah-jumlah itu dan nilai yang dihitung. Anda mungkin dapat membalikkan ekspresi ini untuk mendapatkan estimator pencocokan momen untuk dan - Rosenbaum sebenarnya melakukan itu di bagian 3 dalam kasus terpotong - tetapi bukan itu yang saya inginkan di sini. \E(Z+)\Cov(Z+)LμΣ
Dougal