Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 38
Текст из файла (страница 38)
В нем умножениена обратную матрицу (AP̃ AT + I)−1 заменяется делением на скалярнуювеличину α.24811.7 Рекурсия МНК в стандартной ковариационной формеУпражнение 11.2. С применением леммы 11.1 к выражению (11.24)докажите, что в общем случае ненормализованных экспериментальных данных (11.1) при их статистической интерпретации с ковариацией ошибок v,равной R, матрица Калмана K в алгоритме (11.28) определяется выражениемK = P̃ AT (AP̃ AT + R)−1.(11.32)Упражнение 11.3. Статистическая интерпретация алгоритма (11.31)дана в подразд. 11.3 для случая, когда экспериментальные данные (11.1)нормализованы, то есть характеризуются математическими ожиданиями (11.11). Это выражается добавлением«+1» в выражении для α, (11.31),Tпричем эта «+1» есть не что иное как E vv = 1 для ошибки v в (11.30).Докажите, что в более общем случае, когда матрица R в (11.32) являетсядиагональной,R = diag (r1, r2, . .
. , rm ),(11.33)и только в этом случае матричный алгоритм (11.28) с матрицей K из (11.32)эквивалентен m-кратному повторению скалярного алгоритма вида (11.31) сα = aT P̃ a + r, где aT — i-я строка матрицы A, r = ri и z — i-й элементвектора z, (11.1), при i = 1, 2, . . . , m.Таким образом, скаляризованный алгоритм Калмана (11.31) имеет следующее общее представление:α = aT P̃ a+r,K = P̃ a/α,P̂ = P̃ −KaT P̃ ,x̂ = x̃+K(z −aT x̃), (11.34)если наблюдения (11.1) в их статистической интерпретации составлены изm отдельных, независимых друг от друга, скалярных данных вида (11.30),каждое с ковариацией r = ri, i = 1, 2, . .
. , m.Еще раз отметим, что в применении к решению переопределенной системы алгебраических уравнений в (11.34) следует считать r = 1, то естьиспользовать (11.31).Замечание 11.3.Условие (11.33) не является ни в коей мереограничительным для использования (11.34). Используя разложение Холесского без квадратных корней (см. разд. 6), любую R > 0 можно представить в виде R = U DU T или R = LDLT и затем перейти к измерениямz̄ = U −1 z или z̄ = L−1z, чтобы диагонализировать матрицу ковариацийошибок наблюдений.24911 Оценивание по методу наименьших квадратовУпражнение 11.4.
Докажите, что в алгоритме (11.28) с определениемK по выражению (11.32) вычисление P̂ может быть представлено в так называемой симметричной форме Джозефа:P̂ = (I − KA)P̃ (I − KA)T + KRK T ,которую иначе называют стабилизированным алгоритмом Калмана, таккак она предотвращает возможную потерю положительной определенностиматрицы P̂ , присущую стандартному алгоритму (11.28) с P̂ = P̃ − KAP .При скалярной обработке наблюдений в алгоритме (11.34) это выражениедля P̂ , соответственно, заменяется на скаляризованный алгоритм ДжозефаP̂ = (I − KαT )P̃ (I − αK T ) + rKK T .11.8Ковариационный алгоритм Поттера для МНКВместо матриц P̃ и P̂ , по своей природе положительно определенных,далее оперируем с квадратными корнями из них, соответственно, S̃ и Ŝ,отвечающими равенствам S̃ S̃ T и Ŝ Ŝ T .
Перепишем выражение для P̂ в (11.34)в видеŜ Ŝ T = S̃(In − f f T /α)S̃ T , f = S̃ T a, α = r + f T f,где n — размерность векторов x̂, x̃, и потребуем так выбрать число β, чтобыобеспечить справедливость следующего разложения:In − f f T /α = (In − βf f T )(In − βf f T ).Отсюда для β получаем квадратное уравнение и из двух его решений выбираемpβ = (1/α)/(1 + r/α),поскольку выбор знака ” + ” обеспечивает меньший уровень относительныхошибок при этих вычислениях. Обозначимpγ = 1/(1 + r/α),тогда β = γ/α. В итоге вместо (11.34) получаем следующий ряд формул:f = S̃ T a ,Tα = f f +r,pγ = 1/(1 + r/α) , (11.35)K = S̃f /α ,TŜ = S̃ − γKf ,Tx̂ = x̃ + K(z − a x̃) ,25011.9 Полная статистическая интерпретация рекурсии в МНКкоторый и составляет алгоритм Поттера.
Он численно более устойчив, чемстандартный ковариационный алгоритм Калмана (11.34), но ему эквивалентен. В целом, для него характерно следующее.Вычисление Ŝ в (11.35) равносильно счету с двойной точностью матрицыP̂ в (11.34) при использовании той же разрядности чисел в компьютере,или, иначе, равносильная точность вычислений матрицы P̂ может бытьдостигнута значительно быстрее. Для матрицы P̂ теперь отсутствует опасность потери положительной определенности, присущая операции вычитания в (11.34), поскольку здесь вычисляют Ŝ, а P̂ = Ŝ Ŝ T и P̂ > 0 приdet(Ŝ) 6= 0.
Недостатком алгоритма (11.35) является наличие операцииизвлечения квадратного корня, отдельной для каждого скалярного наблюдения z = aT x+v, и возможная потеря специального (желательно, треугольного) вида матрицы Ŝ в общем случае. Действительно, для экономии памятии объема вычислений обычно стремятся иметь матрицы Ŝ и S̃ в виде треугольных матриц (обе — нижние треугольные или обе — верхние треугольные), что соответствует разложениям Холесского: P̂ = Ŝ Ŝ T и P̃ = S̃ S̃ T .Однако, если стартовать с матрицы S̃ треугольного вида, выполняя инициализацию в соответствии с (11.27), то из-за вычитания в (11.35) матрица Ŝ вобщем случае не остается треугольной.
Например, пусть для Ŝ и S̃ выбранаверхняя треугольная форма. Тогда только при a = λ(1, 0, . . . , 0)T , где λ —скаляр, для Ŝ в (11.35) будет сохранена та же верхняя треугольная форма,благодаря чему выполнение этапа экстраполяции, согласно (11.29), сводитсяк простому присваиванию: S̃ := Ŝ. Если же выбранная для S̃ треугольнаяформа будет утрачена, то этап экстраполяции матрицы потребует предварительной триангуляризации матрицы Ŝ, то есть операции S̃ := triang (Ŝ).Триангуляризация triang (·) должна быть выполнена ортогональным преобразованием матрицы (·), например, преобразованиями Хаусхолдера, Гивенсаили же Грама–Шмидта, которые рассматривались выше (разд.
7).11.9Полная статистическая интерпретация рекурсии вМНКПростая статистическая интерпретация МНК-решения, данная в п. 11.3,позволила объяснить в п. 11.4 идею включения априорных данных в процесс решения задачи о наименьших квадратах. Она была «простой» в томсмысле, что предполагала случайную природу лишь для ошибки v в экспериментальных данных z, (11.1), при необходимости находить апостериорную25111 Оценивание по методу наименьших квадратовvxAAx+zРис. 11.4. Представление экспериментальных данных z в виде результата измерениянеизвестного вектора x с помощью матрицы A в присутствии случайных ошибок vоценку x̂ для некоторого (постоянного) неизвестного вектора x с учетомнекоторой имеющейся (априорной) оценки x̃ (см. рис. 11.4).
Для той интерпретации было достаточно использовать первые два момента случайнойошибки v: математическое ожидание E {v} = 0 и ковариацию E vv T = Pv .Ниже, описывая случайные векторы, мы также ограничиваемся первымидвумя моментами распределения вероятностей, то есть фактически принимаем гипотезу о нормальном (гауссовом) распределении. Отличие будет втом, что с этого момента мы вводим статистическое описание не только дляv, но и для оцениваемого вектора x. При этом надо иметь в виду, что такоеописание для x должно всегда рассматриваться как условное распределение. Именно это обстоятельство делает приводимую ниже статистическуюинтерпретацию МНК полной, т.
е. учитывающей рекурсию процесса обновления оценок в два этапа.1. Экстраполяция оценок по времени между моментами измерений; означает обновление не по измерению, а по времени; его иногда называют«темпоральным обновлением»; эти оценки по отношению к следующему этапу 2 действуют как априорные оценки.2. Обновление оценок по измерениям; означает переход от текущих априорных оценок (текущего априорного условного распределения) к следующим апостериорным оценкам (новому апостериорному условномураспределению вероятностей) для оцениваемого вектора x.Предположим, что ошибка наблюдения v в уравнении (11.1) есть случайный гауссов вектор. Мы предполагаем, что это уравнение нормализовано,т.
е. записано по типу уравнения (11.10), поэтому(11.36)E {v} = 0, E vv T = Im,где E {·} — оператор математического ожидания, 0 — нулевой m-вектор,Im — единичная матрица (m × m)-матрица.25211.9 Полная статистическая интерпретация рекурсии в МНКЗамечание 11.4. Предположение о предварительной нормализациине нарушает общности и введено только для упрощения записей. В любоймомент оно может быть отозвано путем умножения нормализованных данных на корень квадратный L из матрицы R, где LLT = R (см. п. 11.3) иR — положительно определенная матрица (R > 0). Последнее, т.
е. R > 0,означает, что не найдется такого невырожденного преобразования векторанаблюдений z, (11.1), после которого результат содержал бы элементы, свободные от случайных ошибок. Иными словами, все экспериментальные данные содержат случайные погрешности. Случай, когда некоторые измеренияявляются точными, может быть рассмотрен особо [115].Таким образом, плотность распределения вероятностей вектора v определена выражением1−1/2(11.37)fv (ρ) = (2π)m|R |exp − ρT R−1ρ ,2где согласно (11.36), R = E vv T = Im.Допустим вместе с этим, что еще до проведения эксперимента по схемеz = Ax + v(11.38)и независимо от него вы располагаете априорной информацией о неизвестном векторе x. Это предварительное знание, явившееся результатом вашихпредыдущих экспериментов либо просто кем-то другим полученное и вамсообщенное, выражается величиной x̃.
Сама по себе величина x̃ тоже можетбыть не вполне надежной, то есть вам следует ее рассматривать как случайную, а то конкретное ее значение, которое вам сообщено, вы обознача˜ Это ξ˜ имеет смысл реализованного значения случайного вектора x̃.ете ξ.˜ которое вы принимаете за центр распределения вектора x̃,Кроме этого ξ,вам нужно иметь сведения о разбросе величины x̃ относительно центра, т.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
