File Coverage

blib/lib/Statistics/Sequences/Turns.pm

Criterion	Covered	Total	%
statement	66	71	92.9
branch	26	36	72.2
condition	11	14	78.5
subroutine	16	17	94.1
pod	9	9	100.0
total	128	147	87.0

package Statistics::Sequences::Turns;
use 5.008008;
use strict;
use warnings;
use Carp 'croak';
use Statistics::Sequences 0.14;
use base qw(Statistics::Sequences);
$Statistics::Sequences::Turns::VERSION = '0.12';
use Statistics::Zed 0.10;

=pod

=head1 NAME

Statistics::Sequences::Turns - Kendall's turning-points test - of peaks and troughs in a numerical sequence

=head1 SYNOPSIS

 use strict;
 use Statistics::Sequences::Turns 0.12;
 my $turns = Statistics::Sequences::Turns->new();
 $turns->load([2, 0, 8.5, 5, 3, 5.01, 2, 2, 3]); # numbers; or send as "data => $aref" with each stat call
 my $val = $turns->observed(); # or descriptive methods: expected(), variance(), obsdev() and stdev()
 $val = $turns->z_value(); # # or in list context get back both z- and p-value
 $val = $turns->p_value(); # as above, assume data are loaded
 my $href = $turns->stats_hash(values => {observed => 1, p_value => 1}, ccorr => 1); # incl. any other stat-method
 $turns->dump(values => {observed => 1, expected => 1, p_value => 1}, ccorr => 1, flag => 1, precision_s => 3, precision_p => 7);
 # prints: observed = 11.000, expected = 10.900, p_value = 0.5700167

=head1 DESCRIPTION

Implements Kendall's (1973) "turning point test" of sudden changes as peaks and troughs in the values of a numerical sequence. It is sometimes described as a test of "cyclicity", and often used as a test of randomness. However, it simply counts up the number of local maxima and minima within a sequence, regardless of their spacing and magnitude, and so does not indicate if the changes actually cycle between highs and lows, if they are more or less balanced in magnitude, or if any cycling is periodic. Kendall introduced this as a rough test of ups and downs in a sequence ahead of describing more sensitive tests based on autocorrelation and Fourier analysis.

Specifically, for a sequence of numerical data (interval or ordinal) of size I, a count of turns is incremented if the value on trial I, for all I greater than zero and less than I, is, with respect to its immediate neighbours (the values on I - 1 and I + 1), a greater than both neighbours (a peak) or less than both neighbours (a trough). The difference of this observed number from the mean expected number of turns for a randomly generated sequence, taken as a unit of the standard deviation, gives a I-score for assessing the "randomness" of the sequence, i.e., the absence of a factor systematically affecting the frequency of peaks/troughs, given that, for turns, there is "a fairly rapid tendency of the distribution to normality" (Kendall 1973, p. 24).

=head1 METHODS

=head2 new

 $turns = Statistics::Sequences::Turns->new();

Returns a new Turns object. Expects/accepts no arguments but the classname.

=head2 load

 $turns->load(@data);
 $turns->load(\@data);
 $turns->load('sample1' => \@data); # labelled whatever

Loads data anonymously or by name - see L in the Statistics::Data manpage for details on the various ways data can be loaded and then retrieved (more than shown here). Data must be numerical (ordinal, interval type). All elements must be numerical of the method croaks.

=cut

sub load {
    my $self = shift;
    $self->SUPER::load(@_);
    croak __PACKAGE__, '::load All data must be numerical for turns statistics'
      if !$self->all_numeric( $self->access( index => -1 ) );
    return 1;
}

=head2 add, read, unload

See L for these additional operations on data that have been loaded.

=head2 observed

 $v = $turns->observed(); # use anonymously loaded data
 $v = $turns->observed(index => 1); # ... or give the required "index" for the loaded data
 $v = $turns->observed(label => 'mysequence'); # ... or its "label" value
 $v = $turns->observed(data => \@data); # ... or just give the data now

Returns observed number of turns. This is the number of peaks and troughs, starting the count from index 1 of the sequence (a flat array), checking if both its immediate left/right (or past/future) neighbours are lesser than it (a peak) or greater than it (a trough). Wherever the values in successive indices in the sequence are equal, they are treated as a single observation/datum - so the following:

 0 0 1 1 0 1 1 1 0 1

is counted up for turns as

 0 1 0 1 0 1
   * * * *

This shows four turns - two peaks (0 1 0) and two troughs (1 0 1).

Returns 0 if the given list of is empty, or the number of its elements is less than 3.

=cut

sub observed {
    my $self = shift;
    my $args = ref $_[0] ? shift : {@_};
    my $data = _set_data( $self, $args );
    my $trials  = scalar @{$data};
    return 0 if not $trials or $trials < 3;
    my ($count, $i) = (0);
    for ($i = 1; $i < $trials - 1; $i++) {
        if (  ( $data->[ $i - 1 ] > $data->[$i] )
            && ( $data->[ $i + 1 ] > $data->[$i] ) )
        {    # trough at $i
            $count++;
        }
        elsif (( $data->[ $i - 1 ] < $data->[$i] )
            && ( $data->[ $i + 1 ] < $data->[$i] ) )
        {    # peak at $i
            $count++;
        }
    }
    return $count;
}
*turncount_observed = \&observed;
*tco                = \&observed;

=head2 expected

 $v = $turns->expected(); # use first-loaded data; or specify by "index" or "label", or give it as "data" - see observed()
 $v = $turns->expected(data => \@data); # use these data
 $v = $turns->expected(trials => 10); # don't use actual data; calculate from this number of trials

Returns the expected number of turns, which is set by I the number of trials/observations/sample-size ...:

=for html   E[T] = 2 / 3 (N – 2)

or, equivalently (in some sources),

=for html   E[T] = ( 2N – 4 ) / 3

=cut

sub expected {
    my $self = shift;
    my $args = ref $_[0] ? shift : {@_};
    my $trials =
      defined $args->{'trials'}
      ? $args->{'trials'}
      : scalar( @{ _set_data( $self, $args ) } );
   return 2 / 3 * ( $trials - 2 );
   #return (2 * $trials - 4) / 3;
}
*tce                = \&expected;
*turncount_expected = \&expected;

=head2 variance

 $v = $turns->variance(); # use first-loaded data; or specify by "index" or "label", or give it as "data" - see observed()
 $v = $turns->variance(data => \@data); # use these data
 $v = $turns->variance(trials => number); # don't use actual data; calculate from this number of trials

Returns the expected variance in the number of turns for the given length of data I.

=for html   V[T] = (16N – 29 ) / 90

=cut

sub variance {
    my $self = shift;
    my $args = ref $_[0] ? shift : {@_};
    my $trials =
      defined $args->{'trials'}
      ? $args->{'trials'}
      : scalar( @{ _set_data( $self, $args ) } );
    return ( 16 * $trials - 29 ) / 90;
}
*tcv                = \&variance;
*turncount_variance = \&variance;

=head2 obsdev

 $v = $turns->obsdev(); # use data already loaded - anonymously; or specify its "label" or "index" - see observed()
 $v = $turns->obsdev(data => \@data); # use these data

Returns the observed deviation from expectation for the loaded/given sequence: observed I expected turn-count (I - I). Alias of C is supported.

=cut

sub obsdev {
    return observed(@_) - expected(@_);
}
*observed_deviation = \&obsdev;

=head2 stdev

 $v = $turns->stdev(); # use data already loaded - anonymously; or specify its "label" or "index" - see observed()
 $v = $turns->stdev(data => \@data);

Returns square-root of the variance. Aliases C and C (common in other Statistics modules) are supported.

=cut

sub stdev {
    return sqrt variance(@_);
}
*standard_deviation = \&stdev;
*stddev = \&stdev;

=head2 z_value

 $z = $turns->z_value(ccorr => 1); # use data already loaded - anonymously; or specify its "label" or "index" - see observed()
 $z = $turns->z_value(data => $aref, ccorr => 1);
 ($z, $p) = $turns->z_value(data => $aref, ccorr => 1, tails => 2); # same but wanting an array, get the p-value too

Returns the deviation ratio, or I-score, taking the turncount expected from that observed and dividing by the root variance, by default with a continuity correction in the numerator. Called in list context, returns the I-score with its normal distribution, two-tailed I-value.

The data to test can already have been Led, or sent directly as an aref keyed as B.

Optional named arguments B (1 or 2), B (Boolean for the continuity-correction), B (for the statistic, i.e., I-score) and B (for the I-value).

The method can all be called with "sufficient" data: giving, instead of actual data, the B number of turns, and the number of B, the latter being sufficient to compute the expected number of turns and its variance.

Alias C is supported.

=cut

sub z_value {
    my $self = shift;
    my $args = ref $_[0] ? shift : {@_};
    my $data = _set_data( $self, $args );
    my $trials  = defined $args->{'trials'} ? $args->{'trials'} : scalar @{$data};
    my $zed = Statistics::Zed->new();
    my ( $zval, $pval ) = $zed->zscore(
        observed => defined $args->{'observed'}
        ? $args->{'observed'}
        : $self->observed($args),
        expected => $self->expected( trials => $trials ),
        variance => $self->variance( trials => $trials ),
        ccorr => defined $args->{'ccorr'} ? $args->{'ccorr'} : 1,
        tails => $args->{'tails'} || 2,
        precision_s => $args->{'precision_s'},
        precision_p => $args->{'precision_p'},
    );
    return wantarray ? ( $zval, $pval ) : $zval;
}
*tzs              = \&z_value;
*turncount_zscore = \&z_value;
*zscore          = \&z_value;

=head2 p_value

 $p = $turns->p_value(); # using loaded data and default args
 $p = $turns->p_value(ccorr => 0|1, tails => 1|2); # normal-approximation based on loaded data
 $p = $turns->p_value(data => $aref, ccorr => 1, tails => 2); #  using given data (by-passing load and read)

Returns the normal distribution I-value for the deviation ratio (I-score) of the observed number of turns, 2-tailed and continuity-correct by default (or set B => 1 and B => 0, respectively). Other arguments are as for L.

=cut

sub p_value {
    return ( z_value(@_) )[1];
}
*test       = \&p_value;
*turns_test = \&p_value;
*tnt        = \&p_value;

=head2 stats_hash

 $href = $turns->stats_hash(values => {observed => 1, expected => 1, variance => 1, z_value => 1, p_value => 1}, ccorr => 1);

Returns a hashref for the counts and stats as specified in its "values" argument, and with any options for calculating them. See L for details. If calling via a "turns" object, the option "stat => 'turns'" is not needed (unlike when using the parent "sequences" object).

=head2 dump

 $turns->dump(flag => '1|0', text => '0|1|2');

Print test results to STDOUT. See L in the Statistics::Sequences manpage for details.

=cut

sub dump {
    my $self = shift;
    my $args = ref $_[0] ? $_[0] : {@_};
    $args->{'stat'} = 'turns';
    $self->SUPER::dump($args);
    return $self;
}

sub _set_data
{    # Get data via Statistics::Date
     # Remove equivalent successors: e.g., strip 2nd 2 from (3, 2, 2, 7, 2):
    my $self = shift;
    my $args = ref $_[0] ? $_[0] : {@_};
    my $data = $self->access($args)
      ;    # have been already checked to be numerical if previously load()'ed
    ref $data or croak __PACKAGE__, '::Data for counting up turns are needed';
    my @data_u = ();
    for my $i( 0 .. ( scalar @{$data} - 1 ) ) {
        push @data_u, $data->[$i]
          if not scalar @data_u
          or $data->[$i] != $data_u[-1];
    }
    return \@data_u;
}

__END__