File Coverage

lib/Statistics/Sampler/Multinomial.pm

Criterion	Covered	Total	%
statement	72	81	88.8
branch	12	18	66.6
condition	3	9	33.3
subroutine	15	18	83.3
pod	4	4	100.0
total	106	130	81.5

package Statistics::Sampler::Multinomial;

use 5.014;
use warnings;
use strict;

our $VERSION = '0.5';

use Carp;
use Ref::Util qw /is_arrayref/;
use List::Util qw /min sum/;
use List::MoreUtils qw /first_index/;
use Scalar::Util qw /blessed/;

sub new {
    my ($class, %args) = @_;

    my $data = $args{data};
    croak 'data argument not passed'
      if !defined $data;
    croak 'data argument is not an array ref'
      if !is_arrayref ($data);

    my $first_neg_idx = first_index {$_ < 0} @$data;
    croak "negative values passed in data array"
      if $first_neg_idx >= 0;


    my $self = {
        data => $data,
        data_sum_to_one => $args{data_sum_to_one},
    };

    bless $self, $class;

    my $prng = $args{prng};
    $self->_validate_prng_object ($prng);
    $self->{prng}
      =  $prng
      // "${class}::DefaultPRNG"->new;

    return $self;
}

sub _validate_prng_object {
    my ($self, $prng) = @_;

    #  Math::Random::MT::Auto has boolean op overloading
    #  so make sure we don't trigger it or our tests fail
    #  i.e. don't use "if $prng" directly
    #  (and we waste a random number, but that's less of an issue)
    #return 1 if !defined $prng;

    #  no default yet, so croak if not passed
    croak "prng arg is not defined.  "
        . "Math::Random::MT::Auto is a good default to use"
      if not defined $prng;
    croak 'prng arg is not an object'
      if not blessed $prng;
    croak 'prng arg does not have binomial() method'
      if not $prng->can('binomial');

    return 1;
}


sub _initialise {
    my ($self, %args) = @_;

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

    #  caller has not promised they sum to 1
    if ( $self->{data_sum_to_one}) {
        $self->{sum} = 1;
    }
    else {
        my $sum = sum (@$probs);
        if ($sum != 1) {
            my @scaled_probs = map {$_ / $sum} @$probs;
            $probs = \@scaled_probs;
        }
        $self->{sum} = $sum;
    }

    return;
}

sub get_class_count {
    my $self = shift;
    my $aref = $self->{data};
    return scalar @$aref;
}

sub draw {
    my ($self, $args) = @_;

    #  inefficient, but why would one want a single draw?
    #  here for compatibility with SSM::AliasMethod
    return $self->draw_n_samples (1);
}

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

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

    my $data  = $self->{data}
      // croak 'it appears setup has not been run yet';
    my $K    = scalar @$data;
    my $norm = $self->{sum} // do {$self->_initialise; $self->{sum}};

    my @draws;
    my ($sum_p, $sum_n) = (0, 0);

    foreach my $kk (0..($K-1)) {
        if (!$data->[$kk]) {
            $draws[$kk] = 0;
            next;
        }

        my $res = $prng->binomial (
            #  MRMA does not like p>1
            min (1, $data->[$kk] / ($norm - $sum_p)),
            ($n - $sum_n),
        );
        $draws[$kk] = $res;
        $sum_p += $data->[$kk];
        $sum_n += $res;
    }

    return \@draws;
}

#  Cuckoo package to act as a method wrapper
#  to use the perl PRNG stream by default.
# currently croaks because we have no binomial method
package Statistics::Sampler::Multinomial::DefaultPRNG {
    use Carp;
    sub new {
        croak "No default PRNG yet implemented for Statistics::Sampler::Multinomial.\n"
            . "Try Math::Random::MT::Auto.";
        return bless {}, __PACKAGE__;
    }
    #sub rand {
    #    rand();
    #}
    sub binomial {
        ...
    }
    1;
}


1;
__END__