view release on metacpan or  search on metacpan

lib/Statistics/Sampler/Multinomial/AliasMethod.pm  view on Meta::CPAN

package Statistics::Sampler::Multinomial::AliasMethod;

use 5.014;
use warnings;
use strict;

our $VERSION = '1.02';

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/;

use parent qw /Statistics::Sampler::Multinomial/;

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;

    croak 'prng arg is not an object'
      if not blessed $prng;
    croak 'prng arg does not have rand() method'
      if not $prng->can('rand');

    return 1;
}

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

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

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

    #  algorithm and comments stolen/adapted from
    #  https://hips.seas.harvard.edu/blog/2013/03/03/the-alias-method-efficient-sampling-with-many-discrete-outcomes/

    my (@smaller, @larger);
    my @J = (0) x scalar @$probs;
    my @q = (0) x scalar @$probs;
    my $kk = -1;
    my $K = scalar @$probs;

    foreach my $prob (@$probs){
        $kk++;
        $q[$kk] = $K * $prob;
        if ($q[$kk] < 1.0) {
            push @smaller, $kk
        }
        else {
            push @larger, $kk;
        }
    }
    
    # Loop though and create little binary mixtures that
    # appropriately allocate the larger outcomes over the
    # overall uniform mixture.