Неделя 3. Многомерное нормальное распределение и гауссовские процессы

Author

Carlos Buitrago

1 1. Гауссовский случайный вектор

NoteОпределение

Вектор \(\xi=(\xi_1,\ldots,\xi_n)^\top\) называется гауссовским, если его характеристическая функция имеет вид \[ \varphi_\xi(u)=\mathbb E e^{i u^\top\xi}=\exp\left(iu^\top a-\frac12u^\top\Sigma u\right),\quad u\in\mathbb R^n, \] где \(a\in\mathbb R^n\), а \(\Sigma\) — симметричная неотрицательно определённая матрица. Обозначение: \(\xi\sim\mathcal N(a,\Sigma)\).

При этом \(a=\mathbb E\xi\), \(\Sigma_{ij}=\operatorname{cov}(\xi_i,\xi_j)\). Вырожденные нормальные распределения тоже допускаются: например, константа \(c\) имеет закон \(\mathcal N(c,0)\).

ImportantТри эквивалентные характеристики
  1. \(\xi\) — гауссовский вектор.
  2. Существуют матрица \(A\in\mathbb R^{n\times m}\), вектор \(b\in\mathbb R^n\) и независимые \(Z_j\sim\mathcal N(0,1)\) такие, что \(\xi=AZ+b\).
  3. Для любого \(u\in\mathbb R^n\) величина \(u^\top\xi\) нормальна (включая вырожденный случай).

Почему это полезно? Чтобы доказать гауссовость нового вектора, достаточно представить его как детерминированное линейное преобразование исходного гауссовского вектора.

Независимость: если вектор совместно гауссовский, то его координаты независимы тогда и только тогда, когда попарно некоррелированы (\(\Sigma\) диагональна). Без совместной гауссовости это неверно.

Show Python code 
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
rng = np.random.default_rng(2026)
plt.rcParams.update({'figure.figsize': (9, 4.5), 'axes.grid': True, 'grid.alpha': .25})

1.1 Симуляция 1. Корреляция меняет геометрию нормального облака

Сгенерируем \(Z\sim\mathcal N(0,I_2)\) и положим \(X=AZ\). Тогда \(\operatorname{cov}(X)=AA^\top\).

Show Python code 
z = rng.standard_normal((3000, 2))
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
for ax, rho in zip(axes, [-0.85, 0, 0.85]):
    A = np.array([[1., 0.], [rho, np.sqrt(1-rho**2)]])
    x = z @ A.T
    ax.scatter(x[:, 0], x[:, 1], s=5, alpha=.22)
    ax.set(title=fr'$\rho={rho}$', xlabel='$X_1$', ylabel='$X_2$', xlim=(-4,4), ylim=(-4,4))
plt.tight_layout(); plt.show()

1.2 Важная ловушка: нормальные маргиналы не гарантируют совместную нормальность

Пусть \(Z\sim N(0,1)\), \(S\) независимо от \(Z\) и принимает значения \(\pm1\) с вероятностью \(1/2\). Тогда \(X=Z\) и \(Y=SZ\) по отдельности стандартно нормальны и \(\operatorname{cov}(X,Y)=0\), но \(X+Y\) не нормально: с вероятностью \(1/2\) оно равно нулю.

Show Python code 
z = rng.standard_normal(6000)
s = rng.choice([-1, 1], size=z.size)
x, y = z, s*z
fig, ax = plt.subplots(1, 2, figsize=(11, 4))
ax[0].scatter(x, y, s=3, alpha=.25); ax[0].set(title='Нормальные маргиналы, но не совместно нормальны', xlabel='X', ylabel='Y')
ax[1].hist(x+y, bins=90, density=True); ax[1].set(title='Сумма X+Y имеет атом в нуле', xlabel='X+Y')
plt.tight_layout(); plt.show()
print('Эмпирическая ковариация:', np.cov(x, y)[0, 1])

Эмпирическая ковариация: -0.010901520049899338

2 2. Гауссовский процесс

NoteОпределение

Процесс \((X_t)_{t\in T}\) называется гауссовским, если для любого \(n\ge1\) и любых \(t_1,\ldots,t_n\in T\) вектор \[ (X_{t_1},\ldots,X_{t_n}) \] совместно гауссовский.

Важно: мало проверить, что каждый \(X_t\) по отдельности нормально распределён. Нужно проверить совместную нормальность любого конечного набора времён.

Гауссовский процесс задаётся в смысле конечномерных распределений двумя функциями: \[m(t)=\mathbb EX_t,\qquad R(s,t)=\operatorname{cov}(X_s,X_t).\] Для \(t_1,\ldots,t_n\) имеем \[ (X_{t_1},\ldots,X_{t_n})^\top\sim\mathcal N(m_T,K_T),\quad (K_T)_{ij}=R(t_i,t_j). \]

Траектории и распределения: ковариационная функция определяет конечномерные законы, но сама по себе не означает непрерывность траекторий. Например, для независимых стандартных нормальных \(X_t\) при разных \(t\) ковариация равна \(\mathbf 1_{\{s=t\}}\); на непрерывном множестве индексов это не даёт непрерывных траекторий.

3 3. Ковариационная функция и критерий существования

TipНеобходимое и достаточное условие

Функция \(R:T\times T\to\mathbb R\) может быть ковариацией гауссовского процесса тогда и только тогда, когда она

  1. симметрична: \(R(s,t)=R(t,s)\);
  2. неотрицательно определена: для всех \(n,t_1,\dots,t_n,c_1,\dots,c_n\), \[\sum_{i,j=1}^n c_i c_jR(t_i,t_j)\ge0.\]

Функция среднего \(m:T\to\mathbb R\) может быть произвольной.

Доказательство необходимости: \[ \sum_{i,j}c_ic_jR(t_i,t_j) =\operatorname{Var}\left(\sum_i c_iX_{t_i}\right)\ge0. \] Симметрия следует из симметрии ковариации.

Идея достаточности: каждой конечной совокупности времён сопоставим многомерный нормальный закон с вектором средних \((m(t_i))\) и матрицей \((R(t_i,t_j))\). Положительная полуопределённость гарантирует существование этих законов, а их согласованность позволяет применить теорему Колмогорова о продолжении.

3.1 Быстрые проверки кандидата \(R\)

  • Диагональ: \(R(t,t)=\operatorname{Var}(X_t)\ge0\).
  • Неравенство Коши — Буняковского: \(|R(s,t)|^2\le R(s,s)R(t,t)\).
  • Матрицы Грама: для любых выбранных времён все собственные значения \((R(t_i,t_j))\) неотрицательны.

Отрицательная диагональ или отрицательный определитель \(2\times2\) немедленно исключают функцию. Однако проверка лишь нескольких сеток не является строгим доказательством положительной определённости для всех времён.

4 4. Как симулировать гауссовский процесс на сетке

Для \(t_1,\ldots,t_n\) строим \(K=(R(t_i,t_j))\), затем генерируем \(Z\sim N(0,I_n)\) и \(X=m+LZ\), где \(LL^\top=K\). Если \(K\) вырождена, используем спектральное разложение с обрезкой малых отрицательных ошибок округления. На графике изображается только дискретизация процесса; соединение точек линиями не доказывает непрерывность его траекторий.

Show Python code 
def sample_gp(times, mean, kernel, n_paths=5, seed=1):
    local = np.random.default_rng(seed)
    K = np.array([[kernel(s,t) for t in times] for s in times])
    vals, vecs = np.linalg.eigh((K+K.T)/2)
    if vals.min() < -1e-8:
        raise ValueError('Матрица ковариаций не является PSD')
    L = vecs @ np.diag(np.sqrt(np.clip(vals, 0, None)))
    return np.asarray([mean(t) for t in times])[:, None] + L @ local.standard_normal((len(times), n_paths))

times = np.linspace(0, 5, 140)
examples = [
    ('Броуновское движение', lambda s,t: min(s,t)),
    ('Экспоненциальная ковариация', lambda s,t: np.exp(-abs(s-t))),
    ('Квадратично-экспоненциальная', lambda s,t: np.exp(-0.5*(s-t)**2)),
    ('Косинусная', lambda s,t: np.cos(s-t)),
]
fig, axes = plt.subplots(2, 2, figsize=(13, 8))
for ax, (name, kernel) in zip(axes.flat, examples):
    paths = sample_gp(times, lambda t: 0, kernel, 5)
    ax.plot(times, paths); ax.set(title=name, xlabel='t', ylabel='$X_t$')
plt.tight_layout(); plt.show()

5 5. Линейные преобразования процессов

ImportantОсновной принцип

Если \((X_t)\) — гауссовский процесс, а \(f,g,h\) — детерминированные функции (причём \(g(t)\in T\)), то \[Y_t=f(t)X_{g(t)}+h(t)\] тоже гауссовский процесс. Более общо, любая конечная детерминированная линейная комбинация \(X_{g_j(t)}\) гауссова.

Доказательство: для любых \(t_1,\ldots,t_n\) все исходные величины, встречающиеся в формулах \(Y_{t_i}\), составляют конечномерный гауссовский вектор \(V\). Тогда \((Y_{t_1},\ldots,Y_{t_n})^\top=AV+b\) с детерминированными \(A,b\).

Для \(Y_t=f(t)X_{g(t)}+h(t)\): \[m_Y(t)=f(t)m_X(g(t))+h(t),\] \[R_Y(s,t)=f(s)f(t)R_X(g(s),g(t)).\]

Для \(Y_t=\sum_{j=1}^k a_j(t)X_{g_j(t)}\): \[R_Y(s,t)=\sum_{i,j}a_i(s)a_j(t)R_X(g_i(s),g_j(t)).\]

Нелинейные функции обычно разрушают нормальность: при \(X_t\sim N(0,\sigma_t^2)\) и \(\sigma_t>0\) величина \(X_t^2/\sigma_t^2\sim\chi^2_1\), а не нормальна.

6 6. Задача: какие преобразования гауссовы?

Пусть \((X_t,t\ge0)\) — центрированный гауссовский процесс с ковариацией \(r(s,t)\), причём \(r(t,t)>0\) для всех \(t>0\). Определим:

  1. \(Y_t=X_t^2\);
  2. \(Y_t=tX_t\);
  3. \(Y_t=X_{|\sin t|}-tX_t\);
  4. \(Y_t=X_{t^5+t}+X_{t^3}\).

Найти все гауссовские процессы.

1. Нет. При любом \(t>0\), \(X_t\sim N(0,r(t,t))\) невырожденна, а \(X_t^2\) почти наверное неотрицательна и при этом не является константой. Невырожденная нормальная случайная величина принимает отрицательные значения с положительной вероятностью, следовательно \(X_t^2\) не нормальна.

2. Да. Для любых \(t_1,\ldots,t_n\) вектор \((t_1X_{t_1},\ldots,t_nX_{t_n})\) есть диагональное линейное преобразование гауссовского вектора. Ковариация: \(R_Y(s,t)=st\,r(s,t)\).

3. Да. Времена \(|\sin t|\) и \(t\) неслучайны и неотрицательны. Каждый конечный вектор \(Y\) — линейное преобразование конечного гауссовского вектора из значений \(X\) в этих временах. Ковариация: \[R_Y(s,t)=r(|\sin s|,|\sin t|)-t\,r(|\sin s|,t)-s\,r(s,|\sin t|)+st\,r(s,t).\]

4. Да. Времена \(t^5+t\) и \(t^3\) детерминированны и принадлежат \([0,\infty)\). Конечные линейные комбинации значений \(X\) гауссовы. Ковариация: \[R_Y(s,t)=r(s^5+s,t^5+t)+r(s^5+s,t^3)+r(s^3,t^5+t)+r(s^3,t^3).\]

Ответ: 2, 3, 4.

6.1 Симуляция 2. Почему возведение в квадрат ломает нормальность

Для фиксированного \(t>0\) нормируем \(X_t\) так, что \(X_t/\sqrt{r(t,t)}\sim N(0,1)\). Сравним распределения \(Z\), \(Z^2\) и \(tZ\).

Show Python code 
z = rng.standard_normal(100000)
fig, axes = plt.subplots(1, 3, figsize=(14, 3.8))
for ax, (name, data) in zip(axes, [('Z',z), ('Z²',z*z), ('2Z',2*z)]):
    ax.hist(data, bins=90, density=True, alpha=.75)
    ax.set(title=name, xlabel='значение', ylabel='плотность')
plt.tight_layout(); plt.show()

7 7. Задача: какие функции могут быть ковариациями?

Существует ли центрированный гауссовский процесс с указанной ковариацией \(R(s,t)\) на \([0,\infty)\)?

  1. \(R(s,t)=|s-t|\);
  2. \(R(s,t)=\cos(s-t)\);
  3. \(R(s,t)=-st\);
  4. \(R(s,t)=2\).

1. Нет. Для \(s\ne t\) имеем \(R(s,s)=R(t,t)=0\), но \(R(s,t)=|s-t|>0\). Это противоречит \(|R(s,t)|^2\le R(s,s)R(t,t)=0\). Эквивалентно, матрица \(\begin{pmatrix}0&|s-t|\\|s-t|&0\end{pmatrix}\) имеет отрицательное собственное значение.

2. Да. Возьмём независимые \(Z_1,Z_2\sim N(0,1)\) и определим \(X_t=Z_1\cos t+Z_2\sin t\). Тогда \(\mathbb EX_t=0\) и \[\operatorname{cov}(X_s,X_t)=\cos s\cos t+\sin s\sin t=\cos(s-t).\] Вектор любых значений процесса — линейное преобразование \((Z_1,Z_2)\), следовательно процесс гауссов.

3. Нет. Для \(t>0\), \(R(t,t)=-t^2<0\), что невозможно для дисперсии.

4. Да. Возьмём \(Z\sim N(0,2)\) и положим \(X_t=Z\) для всех \(t\). Тогда \(\operatorname{cov}(X_s,X_t)=2\); все конечномерные распределения совместно нормальны, хотя вырождены.

Ответ: 2, 4.

7.1 Симуляция 4. Постоянная ковариация \(R(s,t)=2\)

Если \(R(s,t)=2\), то \(\operatorname{Var}(X_t-X_s)=2+2-2\cdot2=0\), то есть \(X_s=X_t\) почти наверное для каждой фиксированной пары \((s,t)\). На конечной сетке траектории горизонтальны. Утверждение об одновременном равенстве для всех несчётных \(t\) требует отдельно выбирать подходящую версию процесса; построение \(X_t=Z\) такую версию даёт.

Show Python code 
t = np.linspace(0, 5, 100)
fig, ax = plt.subplots(figsize=(9,3.8))
for v in rng.normal(0, np.sqrt(2), 6):
    ax.plot(t, np.full_like(t, v))
ax.set(title='Xₜ = Z,  Z ~ N(0,2)', xlabel='t', ylabel='$X_t$')
plt.show()

8 8. Дополнительные примеры и связь с приложениями

8.1 Броуновское движение

Стандартное броуновское движение — центрированный гауссовский процесс с \(R(s,t)=\min(s,t)\), \(s,t\ge0\), и непрерывными траекториями. Его приращения на непересекающихся интервалах независимы, поскольку совместно гауссовы и некоррелированы.

8.2 Стационарные гауссовские процессы

Для стационарного по распределению центрированного гауссовского процесса ковариация зависит только от разности времён: \(R(s,t)=k(s-t)\). Примеры: \(e^{-|s-t|/\ell}\) и \(\cos(\omega(s-t))\). Не всякая функция разности допустима: нужна положительная определённость.

8.3 Прогнозирование: условное распределение гауссовского вектора

Если \[\begin{pmatrix}X\\Y\end{pmatrix}\sim N\left(\begin{pmatrix}m_X\\m_Y\end{pmatrix},\begin{pmatrix}K_{XX}&K_{XY}\\K_{YX}&K_{YY}\end{pmatrix}\right),\] то при обратимой \(K_{XX}\) \[Y\mid X=x\sim N\left(m_Y+K_{YX}K_{XX}^{-1}(x-m_X),\;K_{YY}-K_{YX}K_{XX}^{-1}K_{XY}\right).\] Именно эта формула лежит в основе Gaussian Process Regression в машинном обучении. При шумных наблюдениях вместо \(K_{XX}\) используют \(K_{XX}+\sigma_\varepsilon^2 I\).

  • ML: байесовская регрессия с гауссовскими процессами, интервалы неопределённости, байесовская оптимизация.
  • Финансы: броуновское движение как модель шума; гауссовские модели процентных ставок и факторов риска.
  • Физика: случайные поля, временные корреляции, тепловые флуктуации.
  • Биология: гладкие латентные траектории в продольных измерениях и временных рядах.

8.4 Симуляция 5. GP-регрессия: прогноз и неопределённость

Используем ядро \(k(s,t)=\exp[-(s-t)^2/(2\ell^2)]\) и шум наблюдения с дисперсией \(\sigma_\varepsilon^2\).

Show Python code 
x_obs = np.linspace(0, 8, 13)
y_obs = np.sin(x_obs) + rng.normal(0, .18, len(x_obs))
x_star = np.linspace(-1, 9, 260)
ell, noise = 1.1, .18
kernel = lambda a,b: np.exp(-.5*((np.asarray(a)[:,None]-np.asarray(b)[None,:])/ell)**2)
K = kernel(x_obs,x_obs) + noise**2*np.eye(len(x_obs))
Ks = kernel(x_obs,x_star)
Kss = kernel(x_star,x_star)
alpha = np.linalg.solve(K, y_obs)
mu = Ks.T @ alpha
V = np.linalg.solve(K, Ks)
cov = Kss - Ks.T @ V
sd = np.sqrt(np.maximum(np.diag(cov),0))
fig, ax = plt.subplots(figsize=(11,4.5))
ax.scatter(x_obs,y_obs,color='black',label='Наблюдения',zorder=5)
ax.plot(x_star,mu,label='Апостериорное среднее')
ax.fill_between(x_star,mu-1.96*sd,mu+1.96*sd,alpha=.25,label='95% интервал для латентной функции')
ax.set(xlabel='t',ylabel='f(t)',title='Gaussian Process Regression'); ax.legend(); plt.show()