Простой способ предсказания следующего значения (марковский предиктор, “муравьиный” алгоритм, рекурсия Левинсона-Дарбина)

#алгоритм

Здравствуйте!
Имеется числовой ряд, программно некоторый массив длины N. Необходим самый простой
алгоритм предсказания следующего значения N+1, основанный на всех предыдущих значениях.
Подскажите пожалуйста такой.
Мой ряд представлен случайными положительными дробными числами. В общем любые несложные
идеи хотелось бы рассмотреть.
UPD::откуда берутся числа?
Есть 50 равноудаленных приемо-передатчиков (ПП), не важно как, но любой может связаться
с любым (полносвязная топология). Есть абстрактные источник и приемник (которые тоже
могут связаться с любым из ПП). Связь между источником и приемником всегда устанавливается
через 5 ПП (источник и приемник естественно не учитываются). В один момент времени,
например раз в сутки, связь (маршрут/путь) переформировывается по некоторым, для нас
неизвестным, случайным (именно случайным, так как, например, некоторые ПП в данный
момент времени недостижимы по какой-то причине) законам. Необходимо предсказать, какой
путь (или 6-10 других вариантов путей) будет выбран следующим основываясь лишь на знании
предыдущих выборов системы.
Много различных алгоритмов я уже перепробовал. Так как нету статистической зависимости
между двумя различными выборами путей (а если и есть, то ее сложно определить, я даже
боюсь утверждать, что смена пути происходит по, например, нормальному распределению
вероятности, но, думаю, что можно предположить это для данной задачи), вот я и пытаюсь
теперь аппроксимировать на основе "весов" маршрутов (произведения "весов" каждой пары
ПП, извлекаемых из матрицы переходов (веса определены на пересечении строки и столбца:
переход из i-го ПП в j-й), построенной по похожему принципу, описанному @northerner,
но более упрощенному).

Так как веса путей вида |1,2,3,4,5 | 1,3,2,4,5 | 5,4,3,2,1|, равны, нехитрые рассчеты
дают чуть больше 2млн. возможных маршрутов. Для каждого из них я вычисляю вес (по матрице).
Аналогично мы можем посчитать вес маршрута для каждого предыдущего уже выбранного системой
маршрута (промежуток - год, для каждого дня в году). Предсказывая поведение значения
веса, мы из двух миллионов выберем 6-10 вариантов наиболее схожих с нашим предсказанным
значением.
На данный момент я дошел лишь до того, что предположив нормальное распределение мы
можем выделить наиболее редко встречающиеся пары в матрице переходов. Например, есть
пары, которых за год еще ни разу не встречалось (у которых в матрице переходов на пересечении
- нуль их приходится заменять так как они портят картину, как я считаю). Тем не менее
эти пары не начинают пока появляться чаще остальных.

Если есть альтернативные идеи/методы предсказания маршрутов, буду рад. Если есть
идеи, как установить зависимость между сменами маршрутов, распределение и другие характеристики,
тоже буду рад.    

Ответы

Ответ 1

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

То есть разбить множество возможных значений [0; 1] на n промежутков (скажем, на
10), подсчитать, сколько значений попадает в каждый и создать случайную величину, у
которой такое же распределение вероятностей.


Ответ 2

Выбор предиктора очень сильно зависит от специфики последовательности. Если количество
разных чисел, которые могут встретиться в последовательности невелико, предлагаю простой
марковский предиктор.
Будем считать, что элементы последовательности принимают значения от 1 до N (иначе
их можно перенумеровать и будет так).
Назовем возможные значения элемнтов состояниями. Пусть мощность множества состояний
(количество разных значений) равна K. Определим матрицу P размером KxK и два одномерных
массива R и C размером K. Получив очередной элемент последовательности X(N), увеличиваем
на единицу значения P(X(N - 1), X(N)), R(X(N - 1)) и C(X(N)). Таким образом постепенно
накапливаются значения: P(I, J) - количество переходов из I в J, R(I) - количество
уходов из I, C(J) - количество приходов в J. P(I, J) / R(I) - оценка вероятности (частота),
находясь в I перейти в J.
Накопив достаточное количество переходов, можно начать "предсказывать": пусть мы
находимся в состоянии Z. Разыгрываем случайную величину, равномерно распределенную
в [0, 1]. Идем по строке Z слева направо, вычитая P(Z, J) / R(Z). Как только получаем
нуль или меньше, стобец, в котором находимся и есть ожидаемый элемент.
Метод хорошо работает при следующих ограничениях:

число состояний действительно невелико;
процесс стационарен (однороден по времени), то есть, если где-то в начале последовательности
из 15 всегда переходили в 42, а в оставшейся части - ни разу, смысла использовать метод нет;
следующий элемент существенно зависит от текущего (и гораздо менее - от совсем старых);
вероятности перехода распределены достаточно неравномерно (в противном случае имеем
белый шум и никакой предиктор ничего не даст).

В описанном алгоритме C(J) оказались не нужны, но иногда требуются, поэтому оставлю и их.


Ответ 3

UPD

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

Модель сигнала

В основе метода - модель авторегрессии (АР) порядка k для выборки i = 0, 1, ..., n-1:  

xi+f0xi-1 + f1xi-2 + f2xi-3 + ... + fk-1xi-k = 0, где i=k, k+1,..., n-1.

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

Тёплицева симметрия

Согласно методу наименьших квадратов, следует минимизировать невязку:
EPS = < (xi + f0xi-1 + f1xi-2 + f2xi-3 + ... + fkxi-k)2 > i = k..n-1,
для чего следует приравнять к нулю частные производные от невязки по коэффициентам
авторегрессии fs для s=0..k-1.
Это приводит к следующей системе уравнений специального вида:

a0f0 + a1f1 + a2f2 + ... + ak-2fk-2 + ak-1fk-1 = -a1,
a1f0 + a0f1 + a1f2 + ... + ak-3fk-2 + ak-2fk-1 = -a2,
a2f0 + a1f1 + a0f2 + ... + ak-4fk-2 + ak-3fk-1 = -a3,
...
ak-2f0 + ak-3f1 + ak-4f2 + ... + a0fk-2 + a1fk-1 = -ak-1,
ak-1f0 + ak-2f1 + ak-3f2 + ... + a1fk-2 + a0fk-1 = -ak,

где a - автокорреляционная функция (АКФ).

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

Вторая особенность - в том, что столбец свободных членов сформирован из  элементов
той же матрицы a, и к такой системе применим алгоритм Левинсона-Дарбина.

Третья особенность - что на самом деле тёплицевой данная матрица является лишь приблизительно.
Например, для модели второго порядка точный вид уравнений такой:

< xi-12 > f0 + < xi-1xi-2 > f1 = < xixi-1 >,
 < xi-1xi-2 > f0 + < xi-22 > f1 = < xixi-2 >,
где индекс i в каждой сумме вида < ui > пробегает значения от 2 до максимума.
Нетрудно заметить, что:


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


Это может привести к неожиданным эффектам для последовательностей с сильным возрастающим
(убывающим) трендом. 

Так, для выборки из первых 10 членов последовательности Фибоначчи
xi = {1, 1, 2, 3, 5, 8, 13, 21, 34, 55}
получается система уравнений вида:
1869 f0 + 1155 f1 = - 3024,
1155 f0 + 714 f1 = - 1869
с правильным решением f0 = -1, f1 = -1, 
но ожидаемой тёплицевой симметрии задачи нет.
При этом использование АКФ даёт систему
4893 f0 + 3024 f1 = - 3024,
3024 f0 + 4893 f1 = - 1869,
использование которой в данном случае неправомерно.
Основной метод борьбы с этими эффектами - центрирование выборки.
Но главное - это проверка возможности серьёзного упрощения модели и применения алгоритма
Дарбина.

Алгоритм Дарбина

Суть идеи Дарбина - поиск решения для матрицы размерности n+1 в виде:
(f'0, f'1, ... , f'n-1, f'n) = (f0, f1, ... , fn-1, 0) + beta (fn-1, fn-2, ... , f0, 1).
Действительно, подстановка этого решения в матрицу размерности n+1 с учётом системы
размерности n даёт:  

-a1 - beta an + anbeta = -a1,
-a2 - beta an-1 + an-1beta = -a2,
...
-an-1 - beta a1 + a1beta = -an,
anf0 + an-1f1 + an-2f2 + ... + a2fn-2 + a1fn-1 + beta(anfn-1 + an-1fn-2 + an-2fn-3
+ ... + a2f1 + a1f0 + a0) = -an+1.
Первые n уравнений системы удовлетворяются при любом значении beta, последнее - при
beta = - (an+1 + anf0 + an-1f1 + an-2f2 + ... + a2fn-2 + a1fn-1)
/ (anfn-1 + an-1fn-2 + an-2fn-3 + ... + a2f1 + a1f0 + a0).  

Функция durbin() возвращает массив векторов авторегрессии для всего диапазона n.
При n=1 f = (-a1/a0), последующие векторы коэффициентов вычисляются рекуррентно (сначала
beta, затем f).

Программная реализация

В программе на языке PHP реализованы следующие функции:


Центрирование массива center().
Скалярное произведение векторов со сдвигом второго вектора scalar_prod().
Вывод массива с текстовым комментарием print_array().
Вывод системы линейных уравнений второго порядка и её решения print_s().
Вычисление автокорреляционной функции (АКФ) acf().
Сравнение точной и тёплицевой систем второго порядка и их решений для заданной выборки
compare_s().
Алгоритм Дарбина для решения системы уравнений специального вида durbin().
Тестирование алгоритма Дарбина test_durbin().


Текст программы

function center(&$arr){
    $len = count($arr);
    $aver = array_sum($arr) / $len; 
    foreach($arr as &$item){
        $item -= $aver;
    }
}

function scalar_prod($a, $b, $shift = 0, &$c = null){
    $scal = 0;
    if(is_null($c)) $cc = []; else $cc = &$c;
    foreach($a as $key => $item){
        $cc[] = $item * $b[$key+$shift];
        $scal += end($cc);
    }
    return  $scal;  
}

function print_array($arr, $str, $n = 11){
    print $str."[";
    foreach($arr as $key => $item){
        if(!(($key+1) % $n)) print "
"; 
        printf ("\"%03d\" => %.3f, ", $key, $item);
    }
    print "]";  
}   

function print_s($a, $b, $str){
    print("

$str");
    printf("
 %.3f f0 + %.3f f1 = %.3f", $a[0][0], $a[0][1], $b[0]);
    printf("
 %.3f f0 + %.3f f1 = %.3f", $a[1][0], $a[1][1], $b[1]);
    $det = $a[0][0]*$a[1][1] - $a[1][0]*$a[0][1];
    $det0 = $b[0]*$a[1][1] - $b[1]*$a[0][1];
    $det1 = $a[0][0]*$b[1] - $a[1][0]*$b[0];
    printf("
Решение: f = [%f, %f]", (float)$det0 / $det, (float)$det1/$det);
}

function acf($flow, $k, $center = -1, $len = null){
    if(is_null($len)){
        $len = count($flow);
    }
    $slice = array_slice($flow, $k, $len-$k);
    for($lag = 0; $lag <= $k; $lag++){
        $result[$lag] = scalar_prod($slice, $flow, $k-$lag);
    }
    if($center != -1){
        $denom = 1.0/$result[0];
        foreach($result as &$res){
            $res *= $denom;
        }
    } 
    return $result;
}

function compare_s($test){
    $m = count($test);  
    $acf2 = acf($test, 0, -1, $m-2);
    $acf1 = acf($test, 1, -1, $m-1);
    $acf = acf($test, 2);
    $a_exact = [ [$acf1[0],$acf1[1]], [$acf1[1],$acf2[0]] ];
    $a = [ [$acf[0],$acf[1]], [$acf[1],$acf[0]] ];
    $b = [-$acf[1], -$acf[2]];
    print_s($a_exact, $b, "Контроль симметрии матрицы

Точная система (порядок
2):");
    print_s($a, $b, "Тёплицева система (порядок 2):");
}

function durbin($acf, $n){
    $ff = [];
    $f = [-$acf[1]/$acf[0]];
    print_array($f, "
f = ");
    $ff[] = $f;
    for($r = 1; $r < $n; $r++){
        $acr = array_reverse(array_slice($acf, 0, $r+1));
        $fr = array_reverse($f);
        $fr[] = 1;
        $f[] = 0;
        $beta = - ($acf[$r+1] + scalar_prod($f, $acr))/scalar_prod($fr, $acr);
        print("
beta[$r] = $beta ");
        $f = array_map(function($a,$b) use($beta){
            return $a+$beta*$b;
        },$f,$fr);
        $ff[] = $f;
        print_array($f, "  f = ");
    }
    return $ff;
}

function test_durbin($arr, $n, $center=0){
    $len = count($arr);     
    if($center){
        center($arr);
    } 
    $eps_arr = 0;
    foreach($arr as $item){
        $eps_arr += $item*$item;
    }
    printf("

Решение по Дарбину (порядок АР = %d, длина выборки = $len, СКО
выборки = %f):", $n, sqrt($eps_arr/($len-1)));
    $a = acf($arr, $n);
    print_array($a, "

АКФ: ");
    $ff = durbin($a,$n);
    $c = array_reverse(end($ff));
    $eps = 0;
    $brr = [];
    for($j=$n; $j<$len; $j++){
        $brr[$j] = $arr[$j]+scalar_prod($c, $arr, $j-$n);
        $eps += pow($brr[$j],2);
    }
    $k = count($f)-1;
    printf("
СКО остатка = %f", sqrt($eps/($len-$n)));
}

$m = 10;
$test = [1.0, 1,0];
for ($i = 2; $i < $m; $i++){
    $test[$i] = $test[$i-1] + $test[$i-2];
}   
print("*** Фибоначчи - 10: ***");
compare_s($test);
test_durbin($test, 2);

print("

*** Фибоначчи-10 с центрированием: ***");
test_durbin($test, 2, 1);

$m=2000;
$dpi = 2*M_PI;
for($j=0; $j<$m; $j++) $arr[$j] = sin($j);// - $dpi*floor($j/$dpi));
print("

*** Авторегрессия по Дарбину (синусная выборка): ***");
compare_s($arr);

test_durbin($arr, 1);
test_durbin($arr, 2);
test_durbin($arr, 3);
test_durbin($arr, 4);
test_durbin($arr, 5);
test_durbin($arr, 6);


Результаты:


*** Фибоначчи - 10: ***

Контроль симметрии матрицы

Точная система (порядок 2):
1869.000 f0 + 1155.000 f1 = -3024.000
1155.000 f0 + 714.000 f1 = -1869.000
Решение: f = [-1.000000, -1.000000]

Тёплицева система (порядок 2):
4893.000 f0 + 3024.000 f1 = -3024.000
3024.000 f0 + 4893.000 f1 = -1869.000
Решение: f = [-0.618007, -0.000030]

Решение по Дарбину (порядок АР = 2, длина выборки = 10, СКО выборки = 23.321426):

АКФ: ["000" => 4893.000, "001" => 3024.000, "002" => 1869.000, ]
f = ["000" => -0.618, ]
beta[1] = -2.98035943134E-5   f = ["000" => -0.618, "001" => -0.000, ]
СКО остатка = 15.285024

*** Фибоначчи-10 с центрированием: ***

Решение по Дарбину (порядок АР = 2, длина выборки = 10, СКО выборки = 17.795443):

АКФ: ["000" => 2496.320, "001" => 1399.520, "002" => 716.420, ]
f = ["000" => -0.561, ]
beta[1] = 0.0398418808262   f = ["000" => -0.583, "001" => 0.040, ]
СКО остатка = 12.412102

*** Авторегрессия по Дарбину (синусная выборка): ***

Контроль симметрии матрицы

Точная система (порядок 2):
999.018 f0 + 539.793 f1 = -539.750
539.793 f0 + 999.016 f1 = 415.712
Решение: f = [-1.080605, 1.000000]

Тёплицева система (порядок 2):
998.969 f0 + 539.750 f1 = -539.750
539.750 f0 + 998.969 f1 = 415.712
Решение: f = [-1.080619, 1.000008]

Решение по Дарбину (порядок АР = 1, длина выборки = 2000, СКО выборки = 0.707169):

АКФ: ["000" => 999.677, "001" => 539.750, ]
f = ["000" => -0.540, ]
СКО остатка = 0.595153

Решение по Дарбину (порядок АР = 2, длина выборки = 2000, СКО выборки = 0.707169):

АКФ: ["000" => 998.969, "001" => 539.750, "002" => -415.712, ]
f = ["000" => -0.540, ]
beta[1] = 1.00000781118   f = ["000" => -1.081, "001" => 1.000, ]
СКО остатка = 0.000009

Решение по Дарбину (порядок АР = 3, длина выборки = 2000, СКО выборки = 0.707169):

АКФ: ["000" => 998.142, "001" => 538.985, "002" => -415.712, "003" => -988.206, ]
f = ["000" => -0.540, ]
beta[1] = 0.999521483275   f = ["000" => -1.080, "001" => 1.000, ]
beta[2] = 0.925724185475   f = ["000" => -0.154, "001" => 0.000, "002" => 0.926, ]
СКО остатка = 0.000488

Решение по Дарбину (порядок АР = 4, длина выборки = 2000, СКО выборки = 0.707169):

АКФ: ["000" => 998.122, "001" => 538.857, "002" => -415.831, "003" => -988.206, "004"
=> -652.029, ]
f = ["000" => -0.540, ]
beta[1] = 0.999342177677   f = ["000" => -1.079, "001" => 0.999, ]
beta[2] = 0.925843078378   f = ["000" => -0.154, "001" => 0.000, "002" => 0.926, ]
beta[3] = 6.00603640925   f = ["000" => 5.406, "001" => 0.000, "002" => 0.000, "003"
=> 6.006, ]
СКО остатка = 0.004358

Решение по Дарбину (порядок АР = 5, длина выборки = 2000, СКО выборки = 0.707169):

АКФ: ["000" => 997.550, "001" => 538.964, "002" => -415.143, "003" => -987.569, "004"
=> -652.029, "005" => 282.984, ]
f = ["000" => -0.540, ]
beta[1] = 0.999977661062   f = ["000" => -1.081, "001" => 1.000, ]
beta[2] = 0.925422594935   f = ["000" => -0.155, "001" => 0.000, "002" => 0.925, ]
beta[3] = 5.96426000454   f = ["000" => 5.364, "001" => 0.000, "002" => 0.000, "003"
=> 5.964, ]
beta[4] = -1.11184318911   f = ["000" => -1.267, "001" => -0.000, "002" => 0.000, "003"
=> 0.000, "004" => -1.112, ]
СКО остатка = 0.000027

Решение по Дарбину (порядок АР = 6, длина выборки = 2000, СКО выборки = 0.707169):

АКФ: ["000" => 996.630, "001" => 538.238, "002" => -415.008, "003" => -986.697, "004"
=> -651.222, "005" => 282.984, "006" => 957.015, ]
f = ["000" => -0.540, ]
beta[1] = 0.999627453765   f = ["000" => -1.080, "001" => 1.000, ]
beta[2] = 0.925654012051   f = ["000" => -0.155, "001" => 0.000, "002" => 0.926, ]
beta[3] = 5.98719502289   f = ["000" => 5.387, "001" => 0.000, "002" => 0.000, "003"
=> 5.987, ]
beta[4] = -1.11131942365   f = ["000" => -1.266, "001" => -0.000, "002" => -0.000, "003"
=> 0.000, "004" => -1.111, ]
beta[5] = -0.877666481891   f = ["000" => -0.291, "001" => -0.000, "002" => -0.000, "003"
=> 0.000, "004" => -0.000, "005" => -0.878, ]
СКО остатка = 0.000360


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


Ответ 4

Ваш массив элементов можно представить как значение некоторой функции в точках 1,
2, 3... N
Тогда, для поиска нового значения (с аргументом функции N+1), можно использовать
интерполирование по Ньютону. Программирует это довольно просто.
Так мы найдём многочлен P(x), который принимает значение a[i] из массива в ячейке
с индексом i (т.е. P(i) = a[i]), тогда, чтобы найти a[N+1] значение, просто подставим
в этот многочлен вместо x, значение N+1. Т.е. ответ = P(N+1).


Ответ 5

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


Ответ 6

Кроме нейронной сети есть еще такой интересный алгоритм: муравьиный :)
Применив к вашей системе, он будет работать как-то так: предположим, что источнику
сигнала при попытке соединения могут дать следующие ответы: подтвеждение готовности
передачи или отказ.
При передаче - увеличиваем "привлекательность" ПП на какую-то величину "А",
при отказе - уменьшаем на "Б".
"А", "Б"- надо найти экспериментальным путем для оптимальной работы системы. Кроме
того надо будет настроить максимумы для счетчиков "привлекательности" и уменьшения
счетчика со временем.