Домашнее задание#

1. Генерация случайных величин с заранее заданной ковариационной матрицей.#

Разложение Холецкого используется для генерации коррелированных между собой случайных величин. Проще говоря, когда есть какой-то набор независимых случайных величин и ковариационная матрица. Как из этого получить набор случайных величин, имеющих такую ковариационную матрицу ?

Это нужно как для моделирования случайных сигналов и физических процессов, так и в качестве вспомогательного элемента других вычислительных методов (Монте-Карло). Решение этой задачи делается с помощью разложения Холецкого. Алгоритм заключается в следующем:

  1. Осуществить разложение Холецкого ковариационной матрицы: \(\boldsymbol{\Sigma}=\mathbf{A} \mathbf{A}^T\)

  2. Сгенерировать случайный вектор \(\mathbf{z}\), компонентами которого являются независимые случайные величины с нормальным распределением

  3. Решением поставленной задачи будет вектор:

\[ \mathbf{x}=\mathbf{m}+\mathbf{A} \mathbf{z} \]

Здесь \(\mathbf{m}\) - это постоянный вектор, составленный из математических ожиданий компонент вектора \(\mathbf{z}\).

Напишите функцию, которая в качестве входного параметра берёт ковариационную матрицу, а возвращает набор случайных величин, действуя по описанному выше алгоритму.

2. Обработка экспериментальных данных..#

  1. Обработайте какую-нибудь лабораторную работу (например, из курса общей физики или просто из папки на диске), требующую проведения прямой по экспериментально полученным точкам. Для решения задачи регрессии рекомендуется использовать библиотеку scikit-learn (sklearn) или scipy.

  2. Создайте прямую с шумом и аналогично обработайте данные.

В обоих пунктах нужно построить график, на который будут нанесены точки и прямая среднеквадратической регрессии.

Данные можно сформировать следующим образом:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from seaborn import load_dataset
from scipy.optimize import curve_fit

# данные из лабораторной работы
data = pd.read_excel("data.xlsx")

T = data["T"][:-1]
U = data["U"][:-1]
I = data["I"][:-1]

# создание прямой с шумом

true_w1 = 2.8
true_w0 = 1.5
xs = np.arange(0, 10, 0.5)
noise = np.random.normal(0, 0.45, size=xs.shape)
ys = true_w1 * xs + true_w0 + noise

3. Правдоподобие для гауссовой вероятностной модели.#

Пусть дана выборка точек на прямой \(\left\{x_i\right\}\).

Максимизируйте правдоподобие (или его логарифм) в гауссовой вероятностной модели:

\[ \prod_i p\left(x_i\right) \rightarrow \max _{\mu, \sigma} \quad p(x)=\frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(x-\mu)^2}{2 \sigma^2}} . \]

4. Гауссовы интегралы для МНК.#

На лекции обсуждался учет влияния систематической погрешности путем усреднения решения задачи МНК по гауссовому нормальному распределению для \(y\)-координат точек выборки: \(\tilde{y}_i \sim \mathcal{N}\left(y_i, s^2\right)\), где погрешность по оси ординат считалась равной \(s\). Обобщите этот вывод на случай, когда каждая точка имеет свою \(y\)-погрешность \(s_i\). Для этого проведите усреднение по многомерному нормальному распределению для \(\tilde{y}_i\) с произвольной симметричной матрицей ковариации \(A^{-1}\) :

\[\begin{split} \begin{array}{r} \tilde{y} \sim \frac{1}{(2 \pi)^{l / 2} \sqrt{\operatorname{det} A^{-1}}} \exp \left(-\frac{(\tilde{y}-y)^T A(\tilde{y}-y)}{2}\right), \\ \text { где } y=\left(\begin{array}{lll} y_i & \ldots & y_l \end{array}\right)^T, \text { a } \tilde{y}=\left(\begin{array}{lll} \tilde{y}_i & \ldots & \tilde{y}_l \end{array}\right)^T . \end{array} \end{split}\]
  1. Покажите, что распределение (2) правильно нормировано. Указание: Выполните замену координат \(\tilde{y}-y=S z\), где матрица \(S\) диагонализует \(A\).

  2. Вычислите неприводимые парные корреляторы \(\left\langle\left\langle\tilde{y}_i \tilde{y}_j\right\rangle\right\rangle\), усредняя по распределению (2). Указание: Сделайте замену \(\tilde{y}-y=Y\). Для вычисления гауссового интеграла с предэкспонентой вычислите интеграл \(\int d^l Y \exp \left(-Y^T A Y / 2+J^T Y\right)\) и выполните дифференцирование по параметрам \(J_i\) (компоненты вектора \(J\) ).

5. Систематические погрешности в МНК.#

Выполните в условиях предыдущей задачи.

  1. Оцените систематические погрешности параметров модели \(w_\alpha\), следуя вычислению, приведенному на лекции, и используя корреляторы, полученные в предыдущем пункте.

  2. Запишите решение в частном случае диагональной матрицы \(A=\operatorname{diag}\left(A_1, \ldots, A_l\right)\). Как следует выбирать величины \(A_i\) для моделирования \(y\)-погрешности \(i\)-ой точки, равной \(s_i\) ?

6*. Систематические погрешности весов в МНК при учёте погрешностей х-координат.#

Выполните оценку погрешности весов \(w_\alpha\), учитывая систематическую погрешность \(x\)-координат точек выборки, усреднив решение задачи МНK по гауссовому нормальному распределению \(\tilde{x}_i \sim \mathcal{N}\left(x_i, s_i^2\right)\). Для простоты считайте погрешности для каждой точки равными: \(s_i=s\) и пренебрегайте погрешностью \(y\)-координат.

Указание: разложите аналитическое решение задачи МНK в ряд Тейлора по отклонениям \(\tilde{x}_i-x_i\), считая такое разложение допустимым.