```
# Copyright 2008, 2009, 2010 Kevin Ryde
# This file is part of Chart.
#
# Chart is free software; you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation; either version 3, or (at your option) any later version.
#
# Chart is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along
# with Chart. If not, see <http://www.gnu.org/licenses/>.
package App::Chart::Series::Calculation;
use 5.010;
use strict;
use warnings;
use Carp;
use List::Util qw(min max);
# use Locale::TextDomain ('App-Chart');
sub identity {
return $_[0];
}
sub delay {
my ($class, $N) = @_;
my @array;
my $pos = $N - 1; # initial extends
return sub {
my ($value) = @_;
my $ret = $array[$pos];
$array[$pos] = $value;
if (++$pos >= $N) { $pos = 0; }
return $ret;
};
}
sub sma_and_stddev {
my ($class, $N) = @_;
my $delay_proc = $class->delay ($N);
my $total = 0;
my $total_squares = 0;
my $count = 0;
return sub {
my ($value) = @_;
# drop old
my $old = $delay_proc->($value);
if (defined $old) {
$total -= $old;
$total_squares -= $old * $old;
} else {
$count++;
}
# add new
$total += $value;
$total_squares += $value * $value;
return ($total / $count,
sqrt(max (0, $total_squares*$count - $total*$total)) / $count);
};
}
sub sum {
my ($class, $N) = @_;
my $delay_proc = $class->delay ($N);
my $total = 0;
return sub {
my ($value) = @_;
# drop old
my $old = $delay_proc->($value);
if (defined $old) {
$total -= $old;
}
# add new
$total += $value;
return $total;
};
}
sub ma_proc_by_weights {
my @weights = @_;
# $a[0] is the newest point, $a[1] the prev, through to $a[$N-1]
my @a;
my $total_weight;
return sub {
my ($value) = @_;
unshift @a, $value; # add new
# keep last $N points
if (@a > @weights) {
pop @a; # drop old
} else {
# new total weight for bigger @a
$total_weight = List::Util::sum (map {$weights[$_]} 0 .. $#a);
}
if ($total_weight == 0) {
return 0;
}
return (List::Util::sum (map {$a[$_] * $weights[$_]} 0 .. $#a)
/ $total_weight);
};
}
#------------------------------------------------------------------------------
# http://mathworld.wolfram.com/LeastSquaresFitting.html
# Least squares generally, including deriving formula using
# derivative==0 as follows:
#
# The sum of squares is
#
# R^2(a,b) = Sum (y[i] - (a + b*x[i]))^2
#
# Partial derivative with b is
#
# d R^2(a,b)
# ---------- = Sum 2 * (y[i]-b*x[i]) * -x[i]
# db
#
# And want that to be zero at the minimum, so
#
# Sum -2*x[i]*y[i] + 2*b*Sum x[i]^2 = 0
#
# Sum x[i]*y[i]
# b = -------------
# Sum x[i]^2
#
# Return a procedure which calculates a linear regression line fit over an
# accumulated window of $N values.
#
# Each call $proc->($y) enters a new value into the window, and the return
# is two values ($a, $b) where the line is $a+$b*X. The last point entered
# is at X=0 and the preceding ones at X=-1, X=-2, etc. A and B are
# #f if not enough data yet.
#
# To prime the window initially, call $proc with $N-1 many points preceding
# the first desired.
#
sub linreg {
my ($class, $N) = @_;
# The X values used are centred around 0,
# $count=4: -1.5, -0.5, 0.5, 1.5
# $count=5: -2, -1, 0, 1, 2
# $count=6: -2.5, -1.5, -0.5, 0.5, 1.5, 2.5
# etc
# But the return is then adjusted $a is based on the last point as X=0
#
# @array,$pos is a circular list of $count many values. The one at
# @array[$pos] is the oldest and is replaced by a new value to cycle that
# in.
#
# $count is how many points are in @array.
#
# $y_total is the total of the past $count many Y values, ie. the values
# in @array.
#
# $xy2_total is the sum of 2*X*Y for each Y value in @array.
#
# $xx2_total is 2 * the sum of X*X for each X value used. This is a
# constant once $count stops at COUNT. A minimum 1 is enforced for the
# degenerate case of $N==0 (there's no slope to in that case but at least
# avoid a divide by zero).
#
my @array;
my $pos = $N - 1; # initial extends
my $count = 0;
my $y_total = 0;
my $xy2_total = 0;
my $xx2_total = 0;
return sub {
my ($y) = @_;
if ($count >= $N) {
# drop oldest point
my $prev = $array[$pos];
$y_total -= $prev;
$xy2_total += ($count-1) * $prev;
} else {
# gaining a point
$count++;
$xy2_total += $y_total; # adj so below is 1 less x
$xx2_total = max (1, linreg_xx2_calc($count));
}
$array[$pos] = $y;
if (++$pos >= $N) { $pos = 0; }
# shift xy products onto 2 less x each
$xy2_total -= ($y_total + $y_total);
# add this point
$y_total += $y;
$xy2_total += ($count-1) * $y;
my $b = $xy2_total / $xx2_total;
return ($y_total/$count + $b*0.5*($count-1),
$b);
};
}
# `xx2-calc' returns 2*Sum(X^2) for the set of N points centred around
# zero as described in linreg-slop-calc-proc below. This means for
# instance,
#
# N=4: 2 * [ (-1.5)^2 + (-0.5)^2 + (0.5)^2 + (1.5)^2 ]
# N=5: 2 * [ (-2)^2 + (-1)^2 + 0^2 + (1)^2 + (2)^2 ]
#
# This is (N^3-N)/6, which can be established by taking successive
# differences or verified by induction (done separately for odds and evens
# is easiest).
#
# N^3-N is always a multiple of 6, since it can be written (N-1)*N*(N+1)
# which is three consecutive numbers so one is certainly a multiple of 3
# and another a multiple of 2. The result is forced to inexact since
# that's what's wanted for the linreg-slope-calc-proc returns.
#
sub linreg_xx2_calc {
my ($N) = @_;
return $N*($N*$N-1) / 6;
}
1;
```