Математические методы

Математическое моделирование физических процессов и математические методы анализа данных являются неотъемлимой частью современной экспериментальной физики. Постоянно возникает потребность как в совершенствовании существующих методов, так и в разработке принципиально новых подходов.

Статистическая регуляризация некорректных обратных задач

Одной из задач, решаемых группой, является популяризация и развитие метода статистической регуляризации, созданного В.Ф. Турчинным в 70-х годах XX века.

Типичной некорректной обратной задачей, возникающей в физике, является уравнение Фредгольма I рода: \[f(y) = \int \limits_a^b dx K(x,y)\varphi(x)\] Фактически, это уравнение описывает следующее: аппаратная функция прибора \(K(x,y)\) действует на иследуемый спектр или иной входной сигнал \(\varphi\), в результате чего исследователь наблюдает выходной сигнал \(f(y)\). Целью исследователя является востановить сигнал \(\varphi\) по известным \(f(y)\) и \(K(x,y)\). Казалось бы, восстановление сигнала не является сложной задачей, поскольку уравнение Фредгольма имеет точное решение. Но уравение Фредгольма некорректно - бесконечно малое изменение начальных условий приводит к конечному изменению решения. Таким образом, наличие шумов, присутсвущих в любом эксперименте, обесценивает попытки решить это уравение точно.

Теория

Расмотрим некую алгебраизацию уравнения Фредгольма: \[f_m = K_{mn}\varphi_n\] С точки зрения математической статистики мы должны должны оценить \(\vec{\varphi}\) по реализации \(\vec{f}\), зная плотность вероятности для \(\vec{f}\) и содержимое матрицы \(K\). Действуя в духе теории принятия решений, мы должны выбрать вектор-функцию \(\vec{S}\), определяющую \(\vec{\varphi}\) на основе \(\vec{f}\) и называемую стратегией. Для того, чтобы определить, какие стратегии более оптимальные, мы введем квадратичную функцию потерь: \[L(\hat{\varphi},\vec{S}) = (\hat{\varphi}-\vec{S})^2,\] где \(\hat{\varphi}\) — наилучшее решение. Согласно баейсовскому подходу рассмотрим \(\vec{\varphi}\) как случайную переменную и переместим нашу неопределенность о \(\vec{\varphi}\) в априорную плотность \(P(\vec{\varphi})\), выражающую достоверность различных возможных законов природы и определяемую на основе информации, сущетсвующей до проведения эксперимента. При таком подходе выбор оптимальной стратегии основывается на минимизации апостериорного риска: \[r_{\vec{S}}(\vec{\varphi}) \equiv E_{\vec{\varphi}}E_{\vec{f}}[L(\vec{\varphi},\vec{S})|\vec{\varphi}]\] Тогда оптимальная стратегия в случае квадратичной функции потерь хорошо изветсна: \[S^{opt} _n= E[\varphi_n|\vec{f}] = \int \varphi_n P(\vec{\varphi}|\vec{f})d\vec{\varphi}\] Апостерионая плотность \(P(\vec{\varphi}|\vec{f})\) определяется по теореме Баейса: \[P(\vec{\varphi}|\vec{f})= \frac{P(\vec{\varphi})P(\vec{f}|\vec{\varphi})}{\int d\vec{\varphi}P(\vec{\varphi})P(\vec{f}|\vec{\varphi})}\] Кроме того, такой подход позволяет определить дисперсию полученного решения: \[\left\langle \sigma_n^2 \right\rangle = \int (\varphi_n - S^{opt}_n)^2 P(\vec{\varphi}|\vec{f})d\vec{\varphi}\] Мы получили решение, введя априорную плотность \(P(\vec{\varphi})\). Можем ли мы сказать, что-либо о том мире функций \(\varphi(x)\), который задается априорной плотностью? Если ответ на этот вопрос отрицательный, то мы должны будем принять все возможные \(\varphi(x)\) равновероятными и вернуться к нерегуляризованному решению. Таким образом, мы должны ответить на этот вопрос положительно. Именно в этом заключается метод статистической регуляризации — регуляризация решения за счет введения дополнительной априорной информации о \(\varphi(x)\). Если исследователь уже обладает какой-либо априорной информацией (априорной плотностью \(P(\vec{\varphi})\)), он может просто вычислить интеграл и получить ответ. В случае, если такой информации нет, в следующем параграфе описывается, какой минимальной информацией может обладать исследователь и как её использовать для получения регулязованного решения.

Априоная информация

Как показали британские ученые, во всем остальном мире любят дифференцировать. Причем, если математик будет задаваться вопросами о правомерности этой операции, то физик оптимистично верит, что законы природы описываются “хорошим” функциями, то есть гладкими. Иначе говоря, он назначает более гладким \(\varphi(x)\) более высокую априорную плотность вероятности. Так давайте попробуем ввести априорную вероятность, основанную на гладкости. Для этого мы вспомним, что введение априорной иформации — это некоторое насилие над миром, принуждающее законы природы выглядеть удобным для нас образом. Это насилие следует свести к минимуму, и, вводя априорную плотность вероятности, необходимо, что бы информация Шеннона относительно \(\varphi(x)\), содержащаяся в \(P(\vec{\varphi})\), была минимальной. Формализуя выше сказанное, выведем вид априорной плотности, основанной на гладкости функции. Для этого мы будем искать условный экстремум информации: \[ I[P(\vec{\varphi})] = \int \ln{P(\vec{\varphi})} P(\vec{\varphi}) d\vec{\varphi} \to min \] При следующих условиях:

  1. Условие на гладкость \(\varphi(x)\). Пусть \(\Omega\) — некоторая матрица, характеризующая гладкость функции. Тогда потребуем, чтобы достигалось определённое значение функционала гладкости: \[\int (\vec{\varphi},\Omega\vec{\varphi}) P(\vec{\varphi}) d\vec{\varphi} = \omega\] Внимательный читатель должен задать вопрос об определении значения параметра \(\omega\). Ответ на него будет дан далее по тексту.

  2. Нормированность вероятности на единицу: \[\int P(\vec{\varphi}) d\vec{\varphi} = 1\]

При этих условиях доставлять минимум функционалу будет следующая функция: \[P_{\alpha}(\vec{\varphi}) = \frac{\alpha^{Rg(\Omega)/2}\det\Omega^{1/2}}{(2\pi)^{N/2}} \exp(-\frac{1}{2} (\vec{\varphi},\alpha\Omega\vec{\varphi}))\] Параметр \(\alpha\) cвязан с \(\omega\), но поскольку у нас нет собственно информации о конкректных значениях функционала гладкости, выяснять, как имеено он связан, бессмысленно. Что же тогда делать с \(\alpha\), спросите вы? Здесь перед вами расрываются три пути: 1. подбирать значение параметра \(\alpha\) вручную, и тем самым перейти к регуляризации Тихонова, 2. усреднить по всем возможным \(\alpha\), предпологая все возможные \(\alpha\) равновероятными, 3. выбрать наиболее вероятное \(\alpha\) по его апостериорной плотности вероятности \(P(\alpha|\vec{f})\). Этот подход верен, если мы прдепологаем, что в экспериментальных данных содержится достаточно ифнормаци об \(\alpha\). Первый случай нам мало интересен. Во втором случае мы получим следующую формулу для решения: \[\left\langle \varphi_i \right\rangle = \frac{\int d\varphi\, \varphi_i P(f|\varphi) \int\limits d\alpha\,P(\alpha) \alpha^{\frac{Rg(\Omega)}{2}} \exp(-\frac{\alpha}{2} (\vec{\varphi},\Omega\vec{\varphi}))}{\int d\varphi P(f|\varphi) \int\limits d\alpha\,P(\alpha) \alpha^{\frac{Rg(\Omega)}{2}} \exp(-\frac{\alpha}{2} (\vec{\varphi},\Omega\vec{\varphi}))}\] Третий случай будет рассмотрен в следующем разделе на примере гауссовых шумов в эксперименте.

Случай гауссовых шумов

Случай, когда ошибки в эксперименте распределны по Гауссу, замечателен тем, что можно получить аналитическое решение нашей задачи. Решение и его ошибка будут иметь следующий вид: \[ \vec{\varphi} = (K^T\Sigma^{-1}K +\alpha^*\Omega)^{-1}K^T\Sigma^{-1^{T}}\vec{f}\] \[ \Sigma_{\vec{\varphi}} = (K^T\Sigma^{-1}K+\alpha^*\Omega)^{-1}\] где \(\Sigma\) - ковариционная матрица многомерного распределения Гаусса, \(\alpha^*\) - наиболее вероятное значение параметра \(\alpha\), которое определяется из условия максимума апостериорной плотности вероятности: \[P(\alpha|\vec{f}) = C' \alpha^{\frac{Rg(\Omega)}{2}}\sqrt{|(K^T\Sigma^{-1}K+\alpha\Omega)^{-1}|}\exp(\frac{1}{2} \vec{f}^T\Sigma^{-1}K^{T}(K^T\Sigma^{-1}K+\alpha\Omega)^{-1}K^T\Sigma^{-1^{T}}\vec{f})\]

В качестве примера рассмотрим востановление спектра, состоящего из двух гауссовых пиков, которые попали под действие интегрального ядра-ступеньки (функции Хевисайда).


Оптимальное планирование эксперимента при помощи функций значимости параметров

Under construction...

Этот раздел дорабатывается...