nabbla (nabbla1) wrote,
nabbla
nabbla1

Big Data, чтоб их...

Решил всё-таки вывести оптимальную обработку N измерений xk, каждое по M компонент (т.е вектор M×1), и на каждое дана своя ковариационная матрица шума измерений Qk размером M×M.

И понял, что лишь для оценки хотя бы первой компоненты нужно взять M×N различных весов. Это хотя бы записать - уже, блин, сложность. В матричном виде - разве что так:



где A1 - матрица весов для оценки первого компонента (всего таких нужно M штук), X - все измерения, собранные воедино.

В моём случае M=6 (три компоненты параллельного переноса и три угла малого поворота - не хочу пока заморачиваться с честными кватернионными вычислениями, для ошибки установки на корабле должно и углов хватить), N=10 для начала.

И выходит: нужно посчитать 6 матриц размером 6×10, или 360 значений!

Одно радует: это не QuatCore должен считать, а персональный компьютер - он это перемелет без каких-либо проблем. Но выкладки какие-то упоротые, и их результат, система из M×(N-1) уравнений с M×(N-1) неизвестными - не лучше...


Исключительно на время выкладок введём разложение Холецкого для матриц Qk:



где L - нижнетреугольная матрица.

Тогда можно выписать, что каждый измеренный вектор xk состоит из искомого значения X и выборки шума:



Здесь ξk - вектора, составленные из независимых случайных величин с нулевым средним и единичной дисперсией. Может, ещё и Гауссовы, если от этого кому-то полегчает. А нужные нам дисперсии и ковариации получаются путём умножения на Lk.

А вместо гигантской матрицы весов A1 мы всё-таки введём N отдельных векторов, каждый размером M×1. Номер компонента X, который мы сейчас хотим оценить (с некими хитрыми весами просуммировав все-все-все результаты измерений, даже по "чужой" компоненте), обзовём индексом z, "чтоб никто не догадался". Оценка будет выглядеть следующим образом:



Это получается просто одно число.

Для того, чтобы оценка получалась несмещённой, должно выполняться условие:



Здесь δ - символ Кронекера, он даёт единицу, если два индекса совпадают и ноль, если не совпадают. В общем, довольно-таки простая идея: по той компоненте, которую пытаешься найти, сумма весов должна быть единичкой, а по остальным - нулевой, чтобы соседние компоненты не "проникли" куда не просят.

Я хотел честно это условие сохранить и решать условную задачу оптимизации, но где-то напортачил, поэтому решил всё-таки упростить себе жизнь - сразу выразить azN через все остальные и решать уже безусловную задачу:



Как-то немного "неаккуратненько", но заметно проще.

Мы уже добились, что математическое ожидание нашей оценки величины равно самой величине, тот самый "несмещённый результат". Вообще, насчёт этой несмещённости много копий сломано. Помню, в "физики шутят" было кривовато переведённое стихотворение про стрелка, который никогда не попадал по мишени, зато стрелял несмещённо - одна пуля левее, другая обязательно правее, одна выше и одна ниже, и нахрена оно нужно - пущай будет смещённое, но близкое к истине! И действительно, много смещённых оценок напридумывали, которые дают дисперсию ниже! Как правило, оценку несколько "зажимают к нулю", когда знают порядок величины. Но это всё кроличья нора, в которую прыгать не очень хочется...

А критерием оптимизации выберем минимум дисперсии результата. Выпишем, чему равна "случайная компонента":


При умножении транспонированного вектора, т.е 1×M на матрицу M×M получаем строчку 1×M, т.е просто каждый компонент берётся с некоторым весом. А благодаря тому что отдельные компоненты ξk независимы, да ещё обладают единичной дисперсией, то дисперсия попросту равна сумме квадратов всех этих весов:







Ну да, вернулись к исходным ковариационным матрицам ВНЕЗАПНО. Осталось это выражение минимизировать, причём варьировать мы можем ВСЕ значения azk, а их тут M×N штук. Но есть M условий, чтобы получить несмещённый результат. Чтобы они "автоматически" выполнились, просто уберём azN:



И теперь просто находим дифференциал от этой штуки, учитывая что это просто число, поэтому можно отдельные слагаемые элементарно транспонировать (как было 1×1, так и останется!), а малые второго порядка выкидываем:





Чтобы получить минимум, необходимо, чтобы все коэффициенты при dazk были нулевыми.

В итоге получаем следующие уравнения:







Как ни странно, они отличаются друг от друга только первым слагаемым.

Но чтобы "по-нормальному" решить систему линейных уравнений, надо все вектора выложить "в один ряд", т.е вместо N-1 векторов M×1 сложить один длиннющий вектор M(N-1)×1. Не знаю, как эту "колбасу" обозначить, пусть bz, потому что Block.

Тогда мы можем наконец-то записать систему уравнений в виде блочной матрицы:



Матрица имеет размер M(N-1)×M(N-1), она умножается на вектор M(N-1)×1, что даёт вектор M(N-1)×1 - его мы составляем из N-1 одинаковых векторов M×1, располагая их "один под другим".

Простой пример

В прошлый раз мы взяли M=2, N=2 и ковариационные матрицы:


(100% корреляция между двумя компонентами, единичная дисперсия),


(единичная дисперсия, коэффициент корреляции, он же в данном случае ковариация: k)

M(N-1) = 2, т.е блочная матрица "вырождается" до одного-единственного верхнего значения, Q1+Q2:



Умножение QN на e1 даст её левый столбец:



И приходим к уравнению:



Да, в прошлый раз, "ручками", мы пришли ровно к этой же системе!

Не знаю, можно ли как-то красиво это всё упростить. Вроде дофига "закономерностей", но как "очевидно" найти всё необходимое - пока не вижу.

А может, стоит здесь разглядеть хорошее итеративное решение, поскольку я ожидаю, что в реальных приложениях очень хорошим "нулевым приближением" будет взять одинаковые веса 1/N по "правильной" координате (ту, которую пытаемся оценить) и нулевые по всем остальным.

Очень уж не хочется "в явном виде" создавать матрицу 54×54, с последующим её преобразованием по методу Гаусса кажется это неправильным.


Теперь бы стряхнуть пыль с компьютерной модели видеоизмерителя параметров сближения, получить с неё десяток измерений с ковариационными матрицами шума измерений - и посмотреть, что же получается. Вот вообще не знаю, чего ожидать!

UPD. Кстати, коэффициенты для оценки всех M компонент выражаются более красиво, чем если брать только одну компоненту под индексом z. Вот такое уравнение:



Слева блочная матрица M(N-1)×M(N-1), набранная из ковариационных матриц шума измерений. Матрица B имеет размеры M(N-1)×M - каждый её столбец i - вереница коэффициентов, нужных для оценки Xi. Справа матрица M(N-1)×M.

Возможно, я капитально туплю, и всё это дело элементарно упрощается...
Tags: Монте-Карло для чайников, математика, работа, странные девайсы
Subscribe

  • А всё-таки есть польза от ковариаций

    Вчера опробовал "сценарий", когда варьируем дальность от 1 метра до 11 метров. Получилось, что грамотное усреднение - это взять с огромными весами…

  • Так есть ли толк в ковариационной матрице?

    Задался этим вопросом применительно к своему прибору чуть более 2 недель назад. Рыл носом землю с попеременным успехом ( раз, два, три, четыре),…

  • Big Data, чтоб их ... (4)

    Наконец-то стряхнул пыль с компьютерной модели сближения, добавил в неё код, чтобы мы могли определить интересующие нас точки, и выписать…

  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 0 comments