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=Φ(ρax−ay1−ρ2−−−−−√)Ryx=Φ(ρay−ax1−ρ2−−−−−√)rxy=1−ρ2−−−−−√2π−−√ϕ(h2−2ρhk+k21−ρ2−−−−−−−−−−−−−√)
Teraz Rosenbaum pokazuje, że:
Pr(V)E[X~∣V]Pr(V)E[X~2∣V]Pr(V)E[X~Y~∣V]=qxRxy+ρqyRyx=Pr(V)+axqxRxy+ρ2ayqyRyx+ρrxy=ρPr(V)+ρaxqxRxy+ρayqyRyx+rxy.(1)(3)(5)
Przydatne będzie również rozważenie specjalnego przypadku (1) i (3) z , tj. Obcinaniem 1d:
ay=−∞
Pr(V)E[X~∣V]Pr(V)E[X~2∣V]=qx=Pr(V)=Qx.(*)(**)
Teraz chcemy rozważyć
[XY]=[μxμy]+[σx00σy][X~Y~]∼N([μXμY],[σ2xρσxσyρσxσyσ2y])=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[X∣X>0]+Pr(X+=0)0=Pr(X>0)(μx+σxE[X~∣X~≥ax])=Qxμx+qxσx,
E[X2+]=Pr(X+>0)E[X2∣X>0]+Pr(X+=0)0=Pr(X~≥ax)E[(μx+σxX~)2∣X~≥ax]=Pr(X~≥ax)E[μ2x+μxσxX~+σ2xX~2∣X~≥ax]=Qxμ2x+qxμxσx+Qxσ2x
Var[X+]=E[X2+]−E[X+]2=Qxμ2x+qxμxσx+Qxσ2x−Q2xμ2x−q2xσ2x−2qxQxμxσx=Qx(1−Qx)μ2x+(1−2Qx)qxμxσx+(Qx−q2x)σ2x.
Aby znaleźć , będziemy potrzebować
Cov(X+,Y+)
E[X+Y+]=Pr(V)E[XY∣V]+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