update

Обновление состояния в реальном времени моделью в пространстве состояний Кальман, фильтрующий

Описание

update эффективно обновляет распределенность в режиме реального времени путем применения одной рекурсии Фильтра Калмана, чтобы вычислить моменты распределенности в течение итогового периода заданных данных об ответе.

Чтобы вычислить моменты распределенности рекурсивным приложением Фильтра Калмана в течение каждого периода в заданных данных об ответе, использовать filter вместо этого.

пример

[nextState,NextStateCov] = update(Mdl,Y) возвращает обновленные моменты распределенности в итоговое время T, обусловленный на распределении текущего состояния, путем применения одной рекурсии Фильтра Калмана к полностью заданной, стандартной модели в пространстве состояний Mdl учитывая T наблюдал ответы Y. nextState и NextStateCov среднее значение и ковариация, соответственно, обновленной распределенности.

пример

[nextState,NextStateCov] = update(Mdl,Y,currentState,CurrentStateCov) инициализирует Фильтр Калмана при распределении текущего состояния со средним currentState и ковариационная матрица CurrentStateCov.

пример

[nextState,NextStateCov] = update(___,Name,Value) дополнительные опции использования, заданные одними или несколькими аргументами name-value и использованием любая из комбинаций входных аргументов в предыдущих синтаксисах. Например, update(Mdl,Y,Params=params,SquareRoot=true) устанавливает неизвестные параметры в частично заданной модели Mdl к значениям в params, и задает использование варианта Фильтра Калмана квадратного корня для числовой устойчивости.

пример

[nextState,NextStateCov,logL] = update(___) также возвращает логарифмическую правдоподобность, вычисленную для каждого наблюдения в Y.

Примеры

свернуть все

Предположим, что скрытый процесс является AR (1). Уравнение состояния

xt=0.5xt-1+ut,

где ut является Гауссовым со средним значением 0 и стандартным отклонением 1.

Сгенерируйте случайную последовательность 100 наблюдений от xt, предположение, что ряд запускается в 1,5.

T = 100;
ARMdl = arima(AR=0.5,Constant=0,Variance=1);
x0 = 1.5;
rng(1); % For reproducibility
x = simulate(ARMdl,T,Y0=x0);

Предположим далее, что скрытый процесс подвергается аддитивной погрешности измерения. Уравнение наблюдения

yt=xt+εt,

где εt является Гауссовым со средним значением 0 и стандартным отклонением 0.75. Вместе, скрытые уравнения процесса и наблюдения составляют модель в пространстве состояний.

Используйте случайный скрытый процесс состояния (x) и уравнение наблюдения, чтобы сгенерировать наблюдения.

y = x + 0.75*randn(T,1);

Задайте четыре содействующих матрицы.

A = 0.5;
B = 1;
C = 1;
D = 0.75;

Задайте модель в пространстве состояний с помощью содействующих матриц.

Mdl = ssm(A,B,C,D)
Mdl = 
State-space model type: ssm

State vector length: 1
Observation vector length: 1
State disturbance vector length: 1
Observation innovation vector length: 1
Sample size supported by model: Unlimited

State variables: x1, x2,...
State disturbances: u1, u2,...
Observation series: y1, y2,...
Observation innovations: e1, e2,...

State equation:
x1(t) = (0.50)x1(t-1) + u1(t)

Observation equation:
y1(t) = x1(t) + (0.75)e1(t)

Initial state distribution:

Initial state means
 x1 
  0 

Initial state covariance matrix
     x1   
 x1  1.33 

State types
     x1     
 Stationary 

Mdl ssm модель. Проверьте, что модель правильно задана с помощью отображения в Командном окне. Программное обеспечение выводит, что процесс состояния является стационарным. Впоследствии, программное обеспечение устанавливает среднее значение начального состояния и ковариацию к среднему значению и отклонению стационарного распределения модели AR (1).

Пропустите наблюдения через модель в пространстве состояний, в режиме реального времени, чтобы получить распределенность в течение периода 100.

[rtfX100,rtfXVar100] = update(Mdl,y)
rtfX100 = 1.2073
rtfXVar100 = 0.3714

update применяет Фильтр Калмана ко всем наблюдениям в y, и возвращает только оценку состояния периода 100.

Сравните результат с результатами filter.

[fX,~,output] = filter(Mdl,y);
size(fX)
ans = 1×2

   100     1

fX100 = fX(100)
fX100 = 1.2073
fXVar100 = output(end).FilteredStatesCov
fXVar100 = 0.3714
tol = 1e-10;
discrepencyMeans = fX100 - rtfX100;
discrepencyVars = fXVar100 - rtfXVar100;
areMeansEqual = norm(discrepencyMeans) < tol
areMeansEqual = logical
   1

areVarsEqual = norm(discrepencyVars) < tol
areVarsEqual = logical
   1

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

Полагайте, что симулированные данные и модель в пространстве состояний в Вычисляют Только Распределение конечного состояния Из Фильтра Калмана.

T = 100;
ARMdl = arima(AR=0.5,Constant=0,Variance=1);
x0 = 1.5;
rng(1); % For reproducibility

x = simulate(ARMdl,T,Y0=x0);

y = x + 0.75*randn(T,1);

A = 0.5;
B = 1;
C = 1;
D = 0.75;
Mdl = ssm(A,B,C,D);

Предположим, что наблюдения доступны последовательно и полагают, что получение обновило обновленную распределенность путем фильтрации каждого нового наблюдения, когда это доступно.

Симулируйте следующую процедуру с помощью цикла.

  1. Создайте переменные, которые хранят моменты распределения начального состояния.

  2. Пропустите входящее наблюдение через модель, задающую текущие моменты распределения начального состояния.

  3. Перезапишите моменты распределения текущего состояния с новыми моментами распределенности.

  4. Повторитесь 2 и 3, когда новые наблюдения доступны.

currentState = Mdl.Mean0;
currentStateCov = Mdl.Cov0;
newState = zeros(T,1);
newStateCov = zeros(T,1);

for j = 1:T
    [newState(j),newStateCov(j)] = update(Mdl,y(j),currentState,currentStateCov);
    currentState = newState(j);
    currentStateCov = newStateCov(j);
end

Постройте наблюдения, истинные значения состояния и новые средние значения состояния каждого периода.

figure
plot(1:T,x,'-k',1:T,y,'*g',1:T,newState,':r','LineWidth',2)
xlabel("Period")
legend(["True state values" "Observations" "New state values"])

Figure contains an axes object. The axes object contains 3 objects of type line. These objects represent True state values, Observations, New state values.

Сравните результаты с результатами filter.

tol = 1e-10;
[fX,~,output] = filter(Mdl,y);
discrepencyMeans = fX - newState;
discrepencyVars = [output.FilteredStatesCov]' - newStateCov;
areMeansEqual = norm(discrepencyMeans) < tol
areMeansEqual = logical
   1

areVarsEqual = norm(discrepencyVars) < tol
areVarsEqual = logical
   1

Фильтр в реальном времени update, примененный целый набор данных последовательно, возвращает те же распределения состояния как filter.

Предположим, что линейное соотношение между изменением в уровне безработицы и темпом роста номинального валового национального продукта (nGNP) представляет интерес. Предположим, что инновациями mismeasured регрессии первого различия уровня безработицы на nGNP темп роста является серия ARMA(1,1) с Гауссовыми воздействиями (то есть, модель регрессии с ARMA (1,1) ошибки и погрешность измерения). Символически, и в форме пространства состояний, модель

[x1,tx2,t]=[ϕθ00][x1,t-1x2,t-1]+[11]utyt-βZt=x1,t+σεt,

где:

  • x1,t ошибочный ряд ARMA в модели регрессии.

  • x2,t фиктивное состояние для MA (1) эффект.

  • y1,t наблюдаемое изменение в уровне безработицы, выкачиваемом темпом роста nGNP (Zt).

  • ut серия Gaussian воздействий, имеющих среднее значение 0 и стандартное отклонение 1.

  • εt серия Gaussian погрешностей измерения со шкалой σ.

Загрузите набор данных Нельсона-Плоссера, который содержит уровень безработицы и nGNP ряд, среди прочего.

load Data_NelsonPlosser

Предварительно обработайте данные путем выполнения этой процедуры:

  1. Удалите продвижение недостающие наблюдения.

  2. Преобразуйте nGNP ряд в ряд возврата при помощи price2ret.

  3. Примените первое различие для ряда уровня безработицы.

vars = ["GNPN" "UR"];
DT = rmmissing(DataTable(:,vars));
T = size(DT,1) - 1; % Sample size after differencing

Z = [ones(T,1) price2ret(DT.GNPN)];
y = diff(DT.UR);

Хотя этот пример удаляет отсутствующие значения, Фильтр Калмана вмещает ряд, содержащий отсутствующие значения.

Задайте содействующие матрицы.

A = [NaN NaN; 0 0];
B = [1; 1];
C = [1 0];
D = NaN;

Задайте модель в пространстве состояний с помощью ssm.

Mdl = ssm(A,B,C,D);

Подбирайте модель ко всем наблюдениям за исключением итоговых 10 наблюдений (хранение выборка). Используйте случайный набор начальных значений параметров для оптимизации. Задайте компонент регрессии и его начальное значение для оптимизации с помощью 'Predictors' и 'Beta0' аргументы пары "имя-значение", соответственно. Ограничьте оценку σ ко всем положительным, вещественным числам.

fh = 10;
params0 = [0.3 0.2 0.2];
[EstMdl,estParams] = estimate(Mdl,y(1:T-fh),params0,'Predictors',Z(1:T-fh,:),...
    'Beta0',[0.1 0.2],'lb',[-Inf,-Inf,0,-Inf,-Inf]);
Method: Maximum likelihood (fmincon)
Sample size: 51
Logarithmic  likelihood:     -87.2409
Akaike   info criterion:      184.482
Bayesian info criterion:      194.141
           |      Coeff       Std Err    t Stat     Prob  
----------------------------------------------------------
 c(1)      |  -0.31780       0.37357    -0.85071  0.39493 
 c(2)      |   1.21242       0.82223     1.47455  0.14034 
 c(3)      |   0.45583       1.32970     0.34281  0.73174 
 y <- z(1) |   1.32407       0.26525     4.99179   0      
 y <- z(2) | -24.48733       1.89161   -12.94520   0      
           |                                              
           |    Final State   Std Dev     t Stat    Prob  
 x(1)      |  -0.38117       0.42842    -0.88971  0.37363 
 x(2)      |   0.23402       0.66222     0.35339  0.72380 

EstMdl ssm модель.

Nowcast уровень безработицы в горизонт прогноза. Симулируйте эту процедуру с помощью цикла:

  1. Вычислите моменты распределения текущего состояния путем пропущения всех наблюдений в выборке через предполагаемую модель.

  2. Когда наблюдение будет доступно в горизонте прогноза, пропустите его через модель. EstMdl не хранит коэффициенты регрессии, таким образом, необходимо передать в них в использовании аргумента Beta значения имени.

  3. Установите моменты состояния распределения текущего состояния на nowcasts.

  4. Повторитесь 2 и 3, когда новые наблюдения будут доступны.

[currentState,currentStateCov] = update(EstMdl,y(1:T-fh),...
    Predictors=Z(1:T-fh,:),Beta=estParams(end-1:end));
unrateF = zeros(fh,2);
unrateCovF = cell(fh,1);

for j = 1:fh
    [unrateF(j,:),unrateCovF{j}] = update(EstMdl,y(T-fh+j),currentState,currentStateCov,...
        Predictors=Z(T-fh+j,:),Beta=estParams(end-1:end));
    currentState = unrateF(j,:)';
    currentStateCov = unrateCovF{j};
end

Постройте предполагаемые, отфильтрованные состояния. Вспомните, что первое состояние является изменением в уровне безработицы, и второе состояние помогает создать первое.

figure
plot(dates((end-fh+1):end),[unrateF(:,1) y((end-fh+1):end)]);
xlabel('Period')
ylabel('Change in the unemployment rate')
title('Filtered Change in the Unemployment Rate')

Figure contains an axes object. The axes object with title Filtered Change in the Unemployment Rate contains 2 objects of type line.

filter функция возвращает только сумму логарифмической правдоподобности для заданных наблюдений. Эффективно вычислить логарифмическую правдоподобность для каждого наблюдения, которое может быть удобным для пользовательских методов оценки, update использования вместо этого.

Полагайте, что симулированные данные и модель в пространстве состояний в Вычисляют Только Распределение конечного состояния Из Фильтра Калмана.

T = 100;
ARMdl = arima(AR=0.5,Constant=0,Variance=1);
x0 = 1.5;
rng(1); % For reproducibility

x = simulate(ARMdl,T,Y0=x0);

y = x + 0.75*randn(T,1);

A = 0.5;
B = 1;
C = 1;
D = 0.75;
Mdl = ssm(A,B,C,D);

Оцените функцию правдоподобия для каждого наблюдения.

[~,~,logLj] = update(Mdl,y);

logL 100 1 вектор; logL(j) логарифмическая правдоподобность, оцененная при наблюдении j.

Используйте filter оценивать вероятность для целого набора данных.

[~,logL] = filter(Mdl,y);

logL скаляр, представляющий полную вероятность данных.

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

tol = 1e-10;
discrepency = logL - sum(logLj);
areEqual = discrepency < tol
areEqual = logical
   1

Входные параметры

свернуть все

Стандартная модель в пространстве состояний в виде ssm объект модели возвращен ssm или estimate.

Если Mdl частично задан (то есть, это содержит неизвестные параметры), задайте оценки неизвестных параметров при помощи 'Params' аргумент значения имени. В противном случае, update выдает ошибку.

Mdl не хранит наблюдаемые ответы или данные о предикторе. Снабдите данными везде, где необходимое использование соответствующего входа или аргументов пары "имя-значение".

Наблюдаемые данные об ответе в виде числовой матрицы или вектора ячейки из числовых векторов.

  • Если Mdl независимо от времени относительно уравнения наблюдения, затем Y T-by-n матрица, где каждая строка соответствует периоду, и каждый столбец соответствует конкретному наблюдению в модели. T является объемом выборки, и m является количеством наблюдений на период. Последняя строка Y содержит последние наблюдения.

  • Если Mdl время, варьируясь относительно уравнения наблюдения, затем Y T-by-1 вектор ячейки. Каждый элемент вектора ячейки соответствует периоду и содержит nt - размерный вектор из наблюдений в течение того периода. Соответствующие размерности содействующих матриц в Mdl.C{t} и Mdl.D{t} должно быть сопоставимо с матрицей в Y{t} в течение всех периодов. Последняя ячейка Y содержит последние наблюдения.

NaN элементы указывают на недостающие наблюдения. Для получения дополнительной информации о том, как Фильтр Калмана вмещает недостающие наблюдения, см. Алгоритмы.

Текущее среднее значение распределенности (другими словами, среднее значение во время 1, прежде чем Фильтр Калмана обрабатывает заданные наблюдения Y) в виде m-by-1 числовой вектор. m является количеством состояний.

Типы данных: double

Текущая ковариационная матрица распределенности (другими словами, ковариационная матрица во время 1, прежде чем Фильтр Калмана обрабатывает заданные наблюдения Y) в виде m-by-m симметричная, положительная полуопределенная числовая матрица.

Типы данных: double

Аргументы name-value

Задайте дополнительные разделенные запятой пары Name,Value аргументы. Name имя аргумента и Value соответствующее значение. Name должен появиться в кавычках. Вы можете задать несколько аргументов в виде пар имен и значений в любом порядке, например: Name1, Value1, ..., NameN, ValueN.

Пример: update(Mdl,Y,Params=params,SquareRoot=true) устанавливает неизвестные параметры в частично заданной модели Mdl к значениям в params, и задает использование варианта Фильтра Калмана квадратного корня для числовой устойчивости.

Оценки неизвестных параметров в частично заданной модели в пространстве состояний MdlВ виде числового вектора.

Если Mdl частично задан (содержит неизвестные параметры, заданные NaNs), необходимо задать Params. estimate функция возвращает оценки параметра Mdl в соответствующей форме. Однако можно предоставить пользовательские оценки путем расположения элементов Params можно следующим образом:

  • Если Mdl явным образом созданная модель (Mdl.ParamMap isempty), расположите элементы Params соответствовать хитам постолбцового поиска NaNs в содействующих матрицах модели в пространстве состояний, векторе средних значений начального состояния и ковариационной матрице.

    • Если Mdl независимо от времени, порядком является ABCD, Mean0, и Cov0.

    • Если Mdl время, варьируясь, порядком является A{1} через A{end}, B{1} через B{end}C1 через C{end}, D{1} через D{end}, Mean0, и Cov0.

  • Если Mdl неявно созданная модель (Mdl.ParamMap указатель на функцию), первый входной параметр функции отображения параметра к матрице определяет порядок элементов Params.

Если Mdl полностью задан, update игнорирует Params.

Пример: Считайте модель в пространстве состояний Mdl с A = B = [NaN 0; 0 NaN] , C = [1; 1], D = 0, и средние значения начального состояния 0 с ковариацией eye(2). Mdl частично задан и явным образом создан. Поскольку параметры модели содержат в общей сложности четыре NaNs, Params должен быть 4 1 вектор, где Params(1) оценка A(1,1), Params(2) оценка A(2,2), Params(3) оценка B(1,1), и Params(4) оценка B(2,2).

Типы данных: double

Отметьте для применения одномерной обработки многомерного ряда (также известный как sequential filtering) в виде true или false. Значение true применяет одномерную обработку.

Одномерная обработка может ускорить и улучшить числовую устойчивость Фильтра Калмана. Однако все инновации наблюдения должны быть некоррелироваными. Таким образом, Dt, Dt' должен быть диагональным, где Dt, t = 1..., T, является одним из следующего:

  • Матричный D{t} в изменяющейся во времени модели в пространстве состояний

  • Матричный D в независимой от времени модели в пространстве состояний

Пример: Univariate=true

Типы данных: логический

Отметьте для применения варианта Фильтра Калмана квадратного корня в виде true или false. Значение true применяет фильтр квадратного корня когда update реализует Фильтр Калмана.

Если вы подозреваете, что собственные значения отфильтрованного состояния или предсказанных ковариационных матриц наблюдения близко к нулю, установите SquareRoot=true. Фильтр квадратного корня устойчив к числовым проблемам, являющимся результатом конечного точность вычислений, но требует большего количества вычислительных ресурсов.

Пример: SquareRoot=true

Типы данных: логический

Предскажите порог неопределенности в виде неотрицательного скаляра.

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

Это - лучшая практика, чтобы установить Tolerance к небольшому числу, например, le-15, преодолеть числовые препятствия во время оценки.

Пример: Tolerance=le-15

Типы данных: double

Переменные предикторы в уравнении наблюдения модели в пространстве состояний в виде T-by-d числовая матрица, где d является количеством переменных предикторов. Строка t соответствует наблюдаемым предикторам в период t (Zt). Расширенное уравнение наблюдения

ytZtβ=Cxt+Dut.

Таким образом, update выкачивает наблюдения с помощью компонента регрессии. β является независимым от времени вектором из коэффициентов регрессии, которые программное обеспечение оценивает всеми другими параметрами.

Если существуют наблюдения n на период, то регрессы программного обеспечения весь ряд предиктора на каждое наблюдение.

Если вы задаете Predictors, затем Mdl должно быть независимо от времени. В противном случае программное обеспечение возвращает ошибку.

По умолчанию программное обеспечение исключает компонент регрессии из модели в пространстве состояний.

Типы данных: double

Коэффициенты регрессии, соответствующие переменным предикторам в виде d-by-n числовая матрица. d является количеством переменных предикторов (см. Predictors).

Если Mdl предполагаемая модель в пространстве состояний, задайте предполагаемые коэффициенты регрессии, сохраненные в estParams.

Выходные аргументы

свернуть все

Среднее значение состояния после update применяет Фильтр Калмана, возвращенный как m-by-1 числовой вектор. Элементы соответствуют порядку состояний, заданных в Mdl (любой строками A или, как определено Mdl.ParamMap).

Ковариационная матрица состояния после update применяет Фильтр Калмана, возвращенный как m-by-m числовая матрица. Строки и столбцы соответствуют порядку состояний, заданных в Mdl (любой строками A или, как определено Mdl.ParamMap).

Логарифмическая правдоподобность для каждого наблюдения в Y, возвращенный как T-by-1 числовой вектор.

Больше о

свернуть все

Обновление распределенности в реальном времени

real-time state-distribution update применяет одну рекурсию Фильтра Калмана к стандартной модели в пространстве состояний, учитывая длину ряд ответа T и распределенность во время T - 1, чтобы вычислить распределенность во время T.

Считайте модель в пространстве состояний описанной в компактной форме

[xtyt]=[AtAtCt]xt1+[Bt0BtCtDt][utεt].

Фильтр Калмана продолжает можно следующим образом в течение каждого периода t:

  1. Получите распределения прогноза в течение каждого периода в данных путем рекурсивного применения условного ожидания к уравнению пространства состояний, учитывая моменты распределения начального состояния x 0|0 и P 0|0, и все наблюдения до времени t − 1 (Yt−11). Получившееся условное распределение

    [xtyt]|Y1t1~Ν([x^t|t1y^t|t1],[Pt|t1Lt|t1Lt|tlVt|t1]),

    где:

    • x^t|t1=Atx^t1|t1, прогноз состояния на время t.

    • y^t|t1=Ctx^t|t1, предсказанный ответ в течение времени t

    • Pt|t1=AtPt1|t1At+BtBt, ковариация прогноза состояния.

    • Vt|t1=CtPt1|t1Ct+DtDt, предсказанная ковариация ответа.

    • Lt|t1=Pt1|t1Ct, состояние и ответ предсказывают ковариацию.

  2. Отфильтруйте наблюдение t через модель, чтобы получить обновленную распределенность:

    xt|Y1t~Ν(x^t|t,Pt|t),

    где:

    • x^t|t=x^t|t1+Lt|t1Vt|t11(yty^t|t1), средство оценки фильтра состояния.

    • Pt|t=Pt|t1Lt|t1Vt|t11Lt|t1, ковариация состояния фильтрует средство оценки.

Когда x^t1|t1 среднее значение текущего состояния и P, t −1|t−1 является ковариацией текущего состояния, x^t|t новое среднее значение состояния и P, t |t является новой ковариацией состояния.

Алгоритмы

  • Фильтр Калмана хранит недостающие данные, не обновляя отфильтрованное оценочное соответствие состояния недостающим наблюдениям. Другими словами, предположите, что существует недостающее наблюдение в период t. Затем прогноз состояния для периода t на основе предыдущего t – 1 наблюдение и отфильтрованное состояние в течение периода t эквивалентен.

  • Для явным образом заданных моделей в пространстве состояний, update применяет все предикторы к каждому ряду ответа. Однако каждый ряд ответа имеет свой собственный набор коэффициентов регрессии.

  • Для КПД, update делает минимальный контроль ввода.

  • В теории ковариационная матрица состояния должна быть симметричной и положительная полуопределенный. update симметрия сил ковариационной матрицы, прежде чем это применит Фильтр Калмана, но это не проверяет, является ли матрица положительной полуопределенный.

Альтернативная функциональность

Чтобы получить отфильтрованные состояния в течение каждого периода в данных об ответе, вызовите filter функцию вместо этого. В отличие от этого, update, filter выполняет всесторонний контроль ввода.

Ссылки

[1] Дербин Дж. и С. Дж. Купмен. Анализ Временных рядов Методами Пространства состояний. 2-й редактор Оксфорд: Издательство Оксфордского университета, 2012.

Введенный в R2021b