Big Data, чтоб их... (2)
Вчера получил упоротое уравнение, чтобы найти, с какими весами нужно брать результаты измерений, чтобы получить наименее шумную и при этом несмещённую оценку вектора параметров:
B=\left(\begin{array}{c}Q_N\\Q_N\\\dots\\Q_N\end{array}\right))
Слева блочная матрица, составленная из ковариационных матриц шума измерений. Для простейшего примера, где N=M=2 (два параметра, измеряются дважды), всё получилось правильно.
Сегодня решил проверить ещё один "примитивный частный случай", когда M=1, то есть мы измеряем одну-единственную величину, но много раз (N раз). Пока ковырял его - нашёл, где ошибся с множителями Лагранжа, решил наконец через них - вышло существенно компактнее...
Сначала всё-таки хочу продемонстрировать, что результат выше тоже приводит к правильным результатам, по крайней мере, в "примитивном частном случае" M=1. Вместо ковариационных матриц получаем просто N дисперсий - значение шума при каждом опыте. Тогда матрица не блочная, а самая обычная, да ещё и можно её вот так представить:
=)
+\left(\begin{array}{c}D_N\\D_N\\\dots\\D_N\end{array}\right)\left(\begin{array}{cccc}1&1&\dots&1\end{array}\right))
И тогда по формуле Шермана-Мориссона её можно относительно легко обратить...
Получается нечто такое (эту здоровенную матрицу назвал C, не знаю почему, надо было её как-то назвать):
-\frac{diag(\frac{1}{D_1},\frac{1}{D_2},\dots,\frac{1}{D_{N-1}})\left(\begin{array}{c}D_N\\D_N\\\dots\\D_N\end{array}\right)\left(\begin{array}{cccc}1&1&\dots&1\end{array}\right)diag(\frac{1}{D_1},\frac{1}{D_2},\dots,\frac{1}{D_{N-1}})}{1+\left(\begin{array}{cccc}1&1&\dots&1\end{array}\right)diag(\frac{1}{D_1},\frac{1}{D_2},\dots,\frac{1}{D_{N-1}})\left(\begin{array}{c}D_N\\D_N\\\dots\\D_N\end{array}\right)})
В знаменателе просто число, и оно оказывается равно
)
И если "увидеть" в единичке дробь DN/DN, получится:
)
Это хороший знак: сумма обратных дисперсий обязана была там вылезти, обзовём её S:
И после некоторых мучений всё-таки приходим к правильному результату:
)
То есть, веса берутся обратно пропорционально дисперсии шума при данном измерении. В сумме же они дают единицу.
Но ведь я тот же результат получал ощутимо проще когда-то!
Записывал дисперсию:
=a_1^2D_1+a_2^2D_2+\dots+a_N^2D_N)
Требовал, чтобы сумма весов была единичной, и находил экстремум с помощью множителей Лагранжа:
=a_1^2D_1+a_2^2D_2+\dots+a_N^2D_N+\lambda(a_1+a_2+\dots+a_N-1))
Занулял все частные производные, что приводило к уравнениям:
Все значения ak выражал через λ подставлял в нижнее уравнение, откуда находил λ, а отсюда и все веса:
В этом "одномерном случае" метод множителей Лагранжа срабатывает на ура!
Кажется, что при замене чисел на вектора и дисперсии на ковариационную матрицу ничего толком не меняется. Но вчера я совершил большую ошибку: нужно было как-то "впихнуть" нулевой вектор в это скалярное выражение (которое пытаемся минимизировать), и я взял квадрат длины этого вектора, после чего получил полнейшую хрень.
Сейчас, смеху ради, попробовал с этим квадратом посмотреть одномерный случай:
=a_1^2D_1+a_2^2D_2+\dots+a_N^2D_N+\lambda(a_1+a_2+\dots+a_N-1)^2)
И тогда частные производные выходят такие:
=0)
Но выражение в скобках заведомо равно нулю, так что λ из всех выкладок "вылетает", и выходят абсурдные результаты
Вся проблема в том, что функция, которая умножается на λ, имеет ВСЕ частные производные равные нулю в интересующих нас точках! По-моему, ни одного такого дурного примера на матане нам не давали, чтобы ответом было "здесь Лагранж неприменим, надо как-то окольными путями идти", так что я и забыл уже о таком ограничении!
Видимо, вместо скалярного λ надо будет взять векторный, т.е первая компонента вектора разности умножить на λ1, вторая компонента на λ2, и так далее:
=\mathbf{a}_{z1}^TQ_1\mathbf{a}_{z1}+\mathbf{a}_{z2}^TQ_2\mathbf{a}_{z2}+\dots+\mathbf{a}_{zN}^TQ_N\mathbf{a}_{zN}+\boldsymbol{\lambda}^T(\mathbf{a}_{z1}+\mathbf{a}_{z2}+\dots+\mathbf{a}_{zN}-\mathbf{e}_z))
Вот теперь, если взять дифференциал и приравнять нулю коэффициенты при всех малых 1-го порядка, получим условия:
И приходим к такому способу решения:
1. Решаем уравнение
\boldsymbol{\lambda}=\mathbf{e}_z)
(z - индекс, показывающий, который из M параметров мы вообще пытаемся найти)
Сумма обратных матриц будет иметь размер M×M, λ - вектор M×1.
2. Находим:
Да, здесь нужно все ковариационные матрицы обратить, но всё-таки обратить десяток матриц 6×6 - это, пожалуй, не так плохо, как решать уравнение с матрицей 54×54...
Можно одним махом находить веса для всех M параметров, что мы пытаемся оценить. Для этого находим матрицу R (от слова Reciprocal - "обратная"):
^{-1})
И затем нужные веса укладываются в N матриц Ak, каждая размером M×M:
Первый столбец Ak указывает, с какими весами нужно брать M параметров, чтобы получить оценку первого параметра, второй столбец - для второго параметра, и т.д.
Такое меня больше устраивает - не нужно эти жуткие блочные матрицы собирать, всё укладывается в размеры M×M, в моём случае 6×6 (три компоненты параллельного переноса и три угла)
Надо, правда, чтобы матрицы были невырождены, так что первый пример со 100% корреляцией здесь не пройдёт. Но "в реальной жизни" у меня 100% корреляции и не будет никогда, так что пущай!
Вот теперь пора это дело испытать на "модельных" результатах...
Слева блочная матрица, составленная из ковариационных матриц шума измерений. Для простейшего примера, где N=M=2 (два параметра, измеряются дважды), всё получилось правильно.
Сегодня решил проверить ещё один "примитивный частный случай", когда M=1, то есть мы измеряем одну-единственную величину, но много раз (N раз). Пока ковырял его - нашёл, где ошибся с множителями Лагранжа, решил наконец через них - вышло существенно компактнее...
Сначала всё-таки хочу продемонстрировать, что результат выше тоже приводит к правильным результатам, по крайней мере, в "примитивном частном случае" M=1. Вместо ковариационных матриц получаем просто N дисперсий - значение шума при каждом опыте. Тогда матрица не блочная, а самая обычная, да ещё и можно её вот так представить:
И тогда по формуле Шермана-Мориссона её можно относительно легко обратить...
Получается нечто такое (эту здоровенную матрицу назвал C, не знаю почему, надо было её как-то назвать):
В знаменателе просто число, и оно оказывается равно
И если "увидеть" в единичке дробь DN/DN, получится:
Это хороший знак: сумма обратных дисперсий обязана была там вылезти, обзовём её S:
И после некоторых мучений всё-таки приходим к правильному результату:
То есть, веса берутся обратно пропорционально дисперсии шума при данном измерении. В сумме же они дают единицу.
Но ведь я тот же результат получал ощутимо проще когда-то!
Записывал дисперсию:
Требовал, чтобы сумма весов была единичной, и находил экстремум с помощью множителей Лагранжа:
Занулял все частные производные, что приводило к уравнениям:
Все значения ak выражал через λ подставлял в нижнее уравнение, откуда находил λ, а отсюда и все веса:
В этом "одномерном случае" метод множителей Лагранжа срабатывает на ура!
Кажется, что при замене чисел на вектора и дисперсии на ковариационную матрицу ничего толком не меняется. Но вчера я совершил большую ошибку: нужно было как-то "впихнуть" нулевой вектор в это скалярное выражение (которое пытаемся минимизировать), и я взял квадрат длины этого вектора, после чего получил полнейшую хрень.
Сейчас, смеху ради, попробовал с этим квадратом посмотреть одномерный случай:
И тогда частные производные выходят такие:
Но выражение в скобках заведомо равно нулю, так что λ из всех выкладок "вылетает", и выходят абсурдные результаты
Вся проблема в том, что функция, которая умножается на λ, имеет ВСЕ частные производные равные нулю в интересующих нас точках! По-моему, ни одного такого дурного примера на матане нам не давали, чтобы ответом было "здесь Лагранж неприменим, надо как-то окольными путями идти", так что я и забыл уже о таком ограничении!
Видимо, вместо скалярного λ надо будет взять векторный, т.е первая компонента вектора разности умножить на λ1, вторая компонента на λ2, и так далее:
Вот теперь, если взять дифференциал и приравнять нулю коэффициенты при всех малых 1-го порядка, получим условия:
И приходим к такому способу решения:
1. Решаем уравнение
(z - индекс, показывающий, который из M параметров мы вообще пытаемся найти)
Сумма обратных матриц будет иметь размер M×M, λ - вектор M×1.
2. Находим:
Да, здесь нужно все ковариационные матрицы обратить, но всё-таки обратить десяток матриц 6×6 - это, пожалуй, не так плохо, как решать уравнение с матрицей 54×54...
Можно одним махом находить веса для всех M параметров, что мы пытаемся оценить. Для этого находим матрицу R (от слова Reciprocal - "обратная"):
И затем нужные веса укладываются в N матриц Ak, каждая размером M×M:
Первый столбец Ak указывает, с какими весами нужно брать M параметров, чтобы получить оценку первого параметра, второй столбец - для второго параметра, и т.д.
Такое меня больше устраивает - не нужно эти жуткие блочные матрицы собирать, всё укладывается в размеры M×M, в моём случае 6×6 (три компоненты параллельного переноса и три угла)
Надо, правда, чтобы матрицы были невырождены, так что первый пример со 100% корреляцией здесь не пройдёт. Но "в реальной жизни" у меня 100% корреляции и не будет никогда, так что пущай!
Вот теперь пора это дело испытать на "модельных" результатах...