Каковы среднее значение и дисперсия многовариантной нормали с 0 цензурой?

9

Пусть находится в . Каковы среднее значение и ковариационная матрица (при этом max вычисляется поэлементно)?ZN(μ,Σ)RdZ+=max(0,Z)

Это происходит, например, потому что, если мы используем функцию активации ReLU внутри глубокой сети и предполагаем через CLT, что входы для данного уровня примерно нормальные, то это распределение выходов.

(Я уверен, что многие люди вычисляли это раньше, но я не мог найти результат, перечисленный где-либо в разумно читаемой форме.)

Дугал
источник
Это упростило бы ваш ответ - возможно, значительно - чтобы заметить, что вы можете получить его, комбинируя результаты двух отдельных вопросов: (1) каковы моменты усеченного нормального распределения и (2) каковы моменты смеси ? Последнее просто и все, что вам нужно сделать, это привести результаты для первого.
whuber
@ Whuber Хмм. Хотя я не сказал этого явно, это, по сути, то, что я делаю в своем ответе, за исключением того, что я не нашел результатов для усеченного двумерного распределения с общим средним значением и дисперсией, и поэтому все равно пришлось выполнять некоторое масштабирование и сдвиг. Есть ли какой-нибудь способ получить, например, ковариацию, не выполняя количество алгебры, которое я должен был сделать? Я, конечно, не утверждаю, что что-то в этом ответе является новым, просто алгебра была утомительной и подверженной ошибкам, и, возможно, кто-то найдет решение полезным.
Дугал
Правильно: я уверен, что ваша алгебра равносильна тому, что я описал, поэтому, похоже, мы разделяем признательность за упрощение алгебры. Один простой способ уменьшить алгебру - это стандартизировать диагональные элементы до единицы, потому что все, что нужно, - это установить единицу измерения для каждой переменной. В этот момент вы можете напрямую включить результаты Розенбаума в (простые, очевидные) выражения для моментов смесей. Может быть, это даже стоит алгебраического упрощения - дело вкуса: без упрощения это приводит к простой модульной компьютерной программе. Σ
whuber
1
Я полагаю, можно написать программу, которая вычисляет моменты непосредственно с результатами Розенбаума и соответствующим образом смешивает их, а затем сдвигает и масштабирует их обратно в исходное пространство. Это, вероятно, было бы быстрее, чем я это сделал.
Дугал

Ответы:

7

Сначала мы можем уменьшить это, чтобы оно зависело только от определенных моментов одномерного / двумерного усеченного нормального распределения: обратите внимание, конечно, что

E[Z+]=[E[(Zi)+]]iCov(Z+)=[Cov((Zi)+,(Zj)+)]ij,
и поскольку мы делаем координатные преобразования некоторых измерений нормального распределения, мы только нужно беспокоиться о среднем значении и дисперсии 1-й цензурированной нормы и ковариации двух 1-мерных норм.

Мы будем использовать некоторые результаты из

С. Розенбаум (1961). Моменты усеченного двумерного нормального распределения . JRSS B, том 23, стр. 405-408. ( Jstor )

Розенбаум считает и рассматривает усечение до события .

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

В частности, мы будем использовать следующие три результата, его (1), (3) и (5). Сначала определите следующее:

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

Теперь Розенбаум показывает, что:

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

Также будет полезно рассмотреть специальный случай (1) и (3) с , то есть 1d усечением: ay=

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

Теперь мы хотим рассмотреть

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

Мы будем использовать которые являются значениями и когда , ,

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

Теперь, используя (*), мы получаем и использование обоих (*) и (**) дает так, чтобы

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.

Чтобы найти , нам понадобится 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,
а затем вычитая мы получаем 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).

Вот код Python для вычисления моментов:

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

и тест Монте-Карло, что он работает:

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

который дает 0.000572145310512 0.00298692620286, указывая, что заявленные ожидания и ковариации соответствуют оценкам Монте-Карло (на основе 10 выборок).10,000,000

Дугал
источник
Вы можете подвести итог, что это за окончательные значения? Это оценки параметров mu и L, которые вы сгенерировали? Может быть, напечатать эти целевые значения?
AdamO
Нет, возвращаемые значения: и ; то, что я напечатал, было расстоянием между оценками Монте-Карло этих величин и вычисленным значением. Вы могли бы, возможно, инвертировать эти выражения, чтобы получить оценку совпадения моментов для и - Розенбаум фактически делает это в своем разделе 3 в усеченном случае - но это не то, что я хотел здесь. \E(Z+)\Cov(Z+)LμΣ
Дугал