В
статье уже выведен алгоритм, но правда линейного предсказания.
Код на С++, но он прокоментирован в статье конечными формулами.
(Оффтоп)
Код:
void ForwardLinearPrediction( vector<double> &coeffs, const vector<double> &x )
{
// GET SIZE FROM INPUT VECTORS
size_t N = x.size() - 1;
size_t m = coeffs.size();
// INITIALIZE R WITH AUTOCORRELATION COEFFICIENTS
vector<double> R( m + 1, 0.0 );
for ( size_t i = 0; i <= m; i++ )
{
for ( size_t j = 0; j <= N - i; j++ )
{
R[ i ] += x[ j ] * x[ j + i ];
}
}
// INITIALIZE Ak
vector<double> Ak( m + 1, 0.0 );
Ak[ 0 ] = 1.0;
// INITIALIZE Ek
double Ek = R[ 0 ];
// LEVINSON-DURBIN RECURSION
for ( size_t k = 0; k < m; k++ )
{
// COMPUTE LAMBDA
double lambda = 0.0;
for ( size_t j = 0; j <= k; j++ )
{
lambda -= Ak[ j ] * R[ k + 1 - j ];
}
lambda /= Ek;
// UPDATE Ak
for ( size_t n = 0; n <= ( k + 1 ) / 2; n++ )
{
double temp = Ak[ k + 1 - n ] + lambda * Ak[ n ];
Ak[ n ] = Ak[ n ] + lambda * Ak[ k + 1 - n ];
Ak[ k + 1 - n ] = temp;
}
// UPDATE Ek
Ek *= 1.0 - lambda * lambda;
}
// ASSIGN COEFFICIENTS
coeffs.assign( ++Ak.begin(), Ak.end() );
}
Хорошо, подсчет автокорреляций явно лишний, вырезаем. Оставляем только рекурсию, в качестве исходных данных подаем предварительно посчитанные ковариации. Пишу в MathCAD, поэтому покорежу С++ код чтобы было похоже.
(Оффтоп)
Код:
void LD(cov)
{
N = length(cov) - 1; // длина вектора исходных данных
m = length(cov); // длина искомых коэффициентов
Ak( m + 1, 0.0 ); // создаем Ak длиной m+1, инициализируем нулями
Ak[ 0 ] = 1.0;
Ek = cov[ 0 ];
// LEVINSON-DURBIN RECURSION
for ( k = 0; k < m; k++ )
{
// COMPUTE LAMBDA
lambda = 0.0;
for ( j = 0; j <= k; j++ )
{
lambda = lambda - Ak[ j ] * R[ k + 1 - j ];
}
lambda = lambda / Ek;
// UPDATE Ak
for ( n = 0; n <= ( k + 1 ) / 2; n++ )
{
temp = Ak[ k + 1 - n ] + lambda * Ak[ n ];
Ak[ n ] = Ak[ n ] + lambda * Ak[ k + 1 - n ];
Ak[ k + 1 - n ] = temp;
}
// UPDATE Ek
Ek = Ek * (1.0 - lambda * lambda);
}
// ASSIGN COEFFICIENTS
Ak; // возвращаем рассчитанные коэффициенты
}
Проверяю встроенной функцией lpcorr расчитывающую PACF.
Код выдает неверные результаты, но если его отлаживать замечаем что при ограничении цикла по k он считает почти верно, но потом затирает найденное...
В приведенной Вами статье я перестаю понимать что происходит на А.4, собственно там теоретический вывод, но где конечные формулы непонятно.
Судя по приведенной мной статье рекурсию можно(?) развернуть в цикл, заведомо создавая векторы максимально необходимой длины, насколько я понял.
Через полные матрицы я научился решать.
Как я уже написал очевидно и в случае Дурбина-Левинсона первые 2 ответа будут:
Напишите пожалуйста как выглядит ответ для следующего x.
это найденный коэффициент в предыдущем шаге.
Хорошо, совпадает.
, где взять 
?
нужны значения с k...