File: RandomWalk.pm

package info (click to toggle)
libstatistics-test-randomwalk-perl 0.02-2
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 124 kB
  • sloc: perl: 199; makefile: 2
file content (296 lines) | stat: -rw-r--r-- 7,298 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
package Statistics::Test::RandomWalk;

use 5.006;
use strict;
use warnings;

our $VERSION = '0.02';

use Carp qw/croak/;
use Params::Util qw/_POSINT _ARRAY _CODE/;
use Memoize;
use Math::BigFloat;
use Statistics::Test::Sequence;
use Class::XSAccessor {
    constructor => 'new',
    getters => {
        rescale_factor => 'rescale',
    },
    setters => {
        set_rescale_factor => 'rescale',
    },
};

=head1 NAME

Statistics::Test::RandomWalk - Random Walk test for random numbers

=head1 SYNOPSIS

  use Statistics::Test::RandomWalk;
  my $tester = Statistics::Test::RandomWalk->new();
  $tester->set_data( [map {rand()} 1..1000000] );
  
  my $no_bins = 10;
  my ($quant, $got, $expected) = $tester->test($no_bins);
  print $tester->data_to_report($quant, $got, $expected);

=head1 DESCRIPTION

This module implements a Random Walk test of a random number generator
as outlined in Blobel et al (Refer to the SEE ALSO section).

Basically, it tests that the numbers C<[0,1]> generated by a
random number generator are distributed evenly. It divides C<[0,1]>
into C<n> evenly sized bins and calculates the number of expected and
actual random numbers in the bin. (In fact, this counts the
cumulated numbers, but that works the same.)

=head1 METHODS

=head2 new

Creates a new random number tester.

=head2 set_rescale_factor

The default range of the random numbers [0, 1) can be rescaled
by a constant factor. This method is the setter for that factor.

=head2 rescale_factor

Returns the current rescaling factor.

=head2 set_data

Sets the random numbers to operate on. First argument
must be either an array reference to an array of random
numbers or a code reference.

If the first argument is a code reference, the second
argument must be an integer C<n>. The code reference is
called C<n>-times and its return values are used as
random numbers.

The code reference semantics are particularly useful if
you do not want to store all random numbers in memory at
the same time. You can write a subroutine that, for example,
generates and returns batches of 100 random numbers so no
more than 101 of these numbers will be in memory at the
same time. Note that if you return 100 numbers at once and
pass in C<n=50>, you will have a sequence of 5000 random
numbers.

=cut

sub set_data {
    my $self = shift;
    my $data = shift;
    if (_ARRAY($data)) {
        $self->{data} = $data;
        return 1;
    }
    elsif (_CODE($data)) {
        $self->{data} = $data;
        my $n = shift;
        if (not _POSINT($n)) {
            croak("'set_data' needs an integer as second argument if the first argument is a code reference.");
        }
        $self->{n} = $n;
        return 1;
    }
    else {
        croak("Invalid arguments to 'set_data'.");
    }
}

=head2 test

Runs the Random Walk test on the data that was previously set using
C<set_data>.

First argument must be the number of bins.

Returns three array references. First is an array of quantiles.
If the number of bins was ten, this (and all other returned arrays)
will hold ten items.

Second are the determined numbers of random numbers below the
quantiles. Third are the expected counts.

=cut

sub test {
    my $self = shift;
    my $bins = shift;
    if (not _POSINT($bins)) {
        croak("Expecting number of bins as argument to 'test'");
    }
    my $rescale_factor = $self->rescale_factor||1;


    my $data = $self->{data};

    if (not defined $data) {
        croak("Set data using 'set_data' first.");
    }

    my $step = 1 / $bins * $rescale_factor;
    my @alpha;
    push @alpha, $_*$step for 1..$bins;

    my @bins = (0) x $bins;
    my $numbers;

    if (_ARRAY($data)) {
        $numbers = @$data;

        foreach my $i (@$data) {
            foreach my $ai (0..$#alpha) {
                if ($i < $alpha[$ai]) {
                    $bins[$_]++ for $ai..$#alpha;
                    last;
                }
            }
        }
    }
    else { # CODE
        my @cache;
        my $calls = $self->{n};
        foreach (1..$calls) {
            # get new data
            push @cache, $data->();
            while (@cache) {
                $numbers++;
                my $this = shift @cache;
                foreach my $ai (0..$#alpha) {
                    if ($this < $alpha[$ai]) {
                        $bins[$_]++ for $ai..$#alpha;
                        last;
                    }
                }
            }
        }
    }

    my @expected_smaller = map Math::BigFloat->new($numbers)*$_/$rescale_factor, @alpha;

    return(
        [map $_/$rescale_factor, @alpha],
        \@bins,
        \@expected_smaller,
    );
}

=head2 data_to_report

From the data returned by the C<test()> method, this
method creates a textual report and returns it as a string.

Do not forget to pass in the data that was returned by C<test()>
or use the C<test_report()> method directly if you do not use
the data otherwise.

=cut

sub data_to_report {
    my $self = shift;
    my $alpha = shift;
    my $got = shift;
    my $expected = shift;
    if (grep {not _ARRAY($_)} ($alpha, $got, $expected)) {
        croak("Please pass the data returned from 'test' to the 'data_to_report' method.");
    }

    my $max_a = _max_length($alpha);
    $max_a = length('Quantile') if length('Quantile') > $max_a;
    my $max_g = _max_length($got);
    $max_g = length('Got') if length('Got') > $max_g;
    my $max_e = _max_length($expected);
    $max_e = length('Expected') if length('Expected') > $max_e;

    my $str = '';

    $str .= sprintf(
        "\%${max_a}s | \%${max_g}s | \%${max_e}s\n",
        qw/Quantile Got Expected/
    );
    $str .= ('-' x (length($str)-1))."\n";
    foreach my $i (0..$#$alpha) {
        $str .= sprintf(
            "\%${max_a}f | \%${max_g}u | \%${max_e}u\n",
            $alpha->[$i], $got->[$i], $expected->[$i]
        );
    }
    return $str;
}

sub _max_length {
    my $max = 0;
    foreach (@{$_[0]}) {
        $max = length $_ if length($_) > $max;
    }
    return $max;

}

=head1 SUBROUTINES

=head2 n_over_k

Computes C<n> over C<k>. Uses Perl's big number support and
returns a L<Math::BigFloat> object.

This sub is memoized.

=cut

memoize('n_over_k');
sub n_over_k {
    my $n = shift;
    my $k = shift;
    my @bits = ((0) x $k, (1) x ($n-$k));
    foreach my $x (1..($n-$k)) {
        $bits[$x-1]--;
    }

    my $o = Math::BigFloat->bone();
    foreach my $i (0..$#bits) {
        $o *= Math::BigFloat->new($i+1)**$bits[$i] if $bits[$i] != 0;
    }

    return $o->ffround(0);
}

1;

__END__

=head1 SEE ALSO

L<Math::BigFloat>, L<Memoize>, L<Params::Util>

Random number generators:
L<Math::Random::MT>, L<Math::Random>, L<Math::Random::OO>,
L<Math::TrulyRandom>, C</dev/random> where available

L<Statistics::Test::Sequence>

The algorithm was taken from: (German)

Blobel, V., and Lohrmann, E. I<Statistische und numerische Methoden
der Datenanalyse>. Stuttgart, Leipzig: Teubner, 1998

=head1 AUTHOR

Steffen Mueller, E<lt>smueller@cpan.orgE<gt>

=head1 COPYRIGHT AND LICENSE

Copyright (C) 2007-2010 by Steffen Mueller

This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself, either Perl version 5.6 or,
at your option, any later version of Perl 5 you may have available.

=cut