File Coverage

blib/lib/Algorithm/Burg.pm

Criterion	Covered	Total	%
statement	59	59	100.0
branch	3	6	50.0
condition	1	6	16.6
subroutine	9	9	100.0
pod	2	2	100.0
total	74	82	90.2

package Algorithm::Burg;
# ABSTRACT: extrapolate time series using Burg's method


use strict;
use warnings qw(all);

use Carp qw(croak);
use List::Util qw(sum);

use Moo;
use MooX::Types::MooseLike::Base qw(
    ArrayRef
    Num
);
use MooX::Types::MooseLike::Numeric qw(
    PositiveInt
);

our $VERSION = '0.001'; # VERSION


has coefficients    => (is => 'rwp', isa => ArrayRef[Num]);


has order           => (is => 'ro', isa => PositiveInt, required => 1);


has series_tail     => (is => 'rwp', isa => ArrayRef[Num]);


sub train {
    my ($self, $time_series) = @_;

    croak '$time_series should be an ArrayRef'
        if ref($time_series) ne 'ARRAY';

    my $m = $self->order;
    my @x = @$time_series;

    croak '$time_series should have more elements than the AR order is'
        if $#x < $m;

    # initialize Ak
    my @Ak = (1.0, (0.0) x $m);

    # initialize f and b
    my @f = @$time_series;
    my @B = @$time_series;

    # Initialize Dk
    my $Dk = sum map {
        2.0 * $f[$_] ** 2
    } 0 .. $#f;
    $Dk -= $f[0] ** 2 + $B[$#x] ** 2;

    # Burg recursion
    for my $k (0 .. $m - 1) {
        # compute mu
        my $mu = sum map {
            $f[$_ + $k + 1] * $B[$_]
        } 0 .. $#x - $k - 1;
        $mu *= -2.0 / $Dk;

        # update Ak
        for my $n (0 .. ($k + 1) / 2) {
            my $t1 = $Ak[$n] + $mu * $Ak[$k + 1 - $n];
            my $t2 = $Ak[$k + 1 - $n] + $mu * $Ak[$n];
            $Ak[$n] = $t1;
            $Ak[$k + 1 - $n] = $t2;
        }

        # update f and b
        for my $n (0 .. $#x - $k - 1) {
            my $t1 = $f[$n + $k + 1] + $mu * $B[$n];
            my $t2 = $B[$n] + $mu * $f[$n + $k + 1];
            $f[$n + $k + 1] = $t1;
            $B[$n] = $t2;
        }

        # update Dk
        $Dk = (1.0 - $mu ** 2) * $Dk
            - $f[$k + 1] ** 2
            - $B[$#x - $k - 1] ** 2;
    }

    $self->_set_series_tail([ @x[$#x - $m .. $#x] ]);
    return $self->_set_coefficients([ @Ak[1 .. $#Ak] ]);
}


sub predict {
    my ($self, $n) = @_;

    my $coeffs = $self->coefficients;
    my $m = $self->order;
    $n ||= $m
        if !$n || $n > $m;

    my @predicted = @{ $self->series_tail };
    for my $i ($m .. $m + $n) {
        $predicted[$i] = -1.0 * sum map {
            $coeffs->[$_] * $predicted[$i - 1 - $_]
        } 0 .. $m - 1;
    }

    return [ @predicted[$m .. $#predicted] ];
}


1;

__END__