File Coverage

blib/lib/Statistics/Distribution/Generator.pm

Criterion	Covered	Total	%
statement	133	199	66.8
branch	17	34	50.0
condition	6	19	31.5
subroutine	38	53	71.7
pod	7	7	100.0
total	201	312	64.4

package Statistics::Distribution::Generator;

use 5.010;

BEGIN {
  if ($] lt '5.012') {
    use strict;
    use warnings;
  }
}

use overload (
  '0+' => '_render',
  '""' => '_render',
  '@{}' => '_render',
  '|' => '_add_alternative',
  'x' => '_add_dimension',
  fallback => 1,
);

use List::AllUtils qw( reduce );
use Exporter qw( import );
use vars qw( $VERSION );

$VERSION = '1.003';

sub logistic ();

our @EXPORT_OK = qw( gaussian uniform logistic supplied gamma exponential dice );
our %EXPORT_TAGS = (all => \@EXPORT_OK);

our $pi = 3.1415926535897932384626433832795028841971693993751;
our $two_pi = 2 * $pi;
our $e = exp 1;

sub _render {
  my $self = shift;
  if ($self->{ dims }) {
    return [ map { $_->_render } @{$self->{ dims }} ];
  }
  elsif ($self->{ alts }) {
    my $accum = reduce { $a + $b } map { $_->{ weight } // 1 } @{$self->{ alts }};
    my $n = rand() * $accum;
    my $answer;
    for my $alt (@{$self->{ alts }}) {
      $n -= ($alt->{ weight } // 1);
      if ($n <= 0) {
        $answer = $alt->_render;
        last;
      }
    }
    return $answer;
  }
  die("Can't render a(n) " . (ref($self)||$self));
}

sub gaussian {
  my ($mean, $sigma) = @_;
  $mean //= 0;
  $sigma //= 1;
  return bless { mean => $mean, sigma => $sigma }, 'Statistics::Distribution::Generator::gaussian';
}

sub uniform {
  my ($min, $max) = @_;
  $min //= 0;
  $max //= 1;
  return bless { min => $min, max => $max }, 'Statistics::Distribution::Generator::uniform';
}

sub logistic () {
  return bless { }, 'Statistics::Distribution::Generator::logistic';
}

sub supplied {
  my ($iv) = @_;
  my $rv;
  if (ref $iv eq 'CODE') {
    $rv = { code => $iv };
  }
  else {
    $rv = { code => sub { return $iv } };
  }
  return bless $rv, 'Statistics::Distribution::Generator::supplied';
}

sub gamma {
  my ($order, $scale) = map { $_ // 1 } @_;
  return bless {
    order => $order,
    scale => $scale,
    norder => int($order),
  }, 'Statistics::Distribution::Generator::gamma';
}

sub exponential {
  my ($lambda) = map { $_ // 1 } @_;
  return bless { lambda => $lambda }, 'Statistics::Distribution::Generator::exponential';
}

sub dice {
  my ($numdice, $numsides) = @_;
  return bless {
    numdice => $numdice,
    numsides => $numsides,
  } => 'Statistics::Distribution::Generator::dice';
}

sub _rand_nonzero {
  my $rv;
  1 while (!($rv = rand));
  return $rv;
}

sub _add_alternative {
  my ($lhs, $rhs, $swapped) = @_;
  ($lhs, $rhs) = ($rhs, $lhs) if $swapped;
  $rhs = supplied($rhs) unless ref($rhs) =~ /^Statistics::Distribution::Generator/;
  if (!exists($lhs->{ alts })) {
    $lhs = {
      alts => [ $lhs ],
    };
  }
  if (!exists($rhs->{ alts })) {
    $rhs = {
      alts => [ $rhs ],
    };
  }
  push @{$lhs->{ alts }}, @{$rhs->{ alts }};
  return bless $lhs, 'Statistics::Distribution::Generator';
}

sub _add_dimension {
  my ($lhs, $rhs, $swapped) = @_;
  ($lhs, $rhs) = ($rhs, $lhs) if $swapped;
  $rhs = supplied($rhs) unless ref($rhs) =~ /^Statistics::Distribution::Generator/;
  if (!exists($lhs->{ dims })) {
    $lhs = {
      dims => [ $lhs ],
    };
  }
  if (!exists($rhs->{ dims })) {
    $rhs = {
      dims => [ $rhs ],
    };
  }
  push @{$lhs->{ dims }}, @{$rhs->{ dims }};
  return bless $lhs, 'Statistics::Distribution::Generator';
}

1;

package Statistics::Distribution::Generator::gaussian;

use 5.010;
use base qw( Statistics::Distribution::Generator );
use overload (
  '0+' => '_render',
  '""' => '_render',
  '|' => '_add_alternative',
  'x' => '_add_dimension',
  fallback => 1,
);

sub _render {
  my $self = shift;
  my $U = rand;
  my $V = rand;
  return $self->{ mean } + (sqrt(-2 * log $U) * cos($two_pi * $V) * $self->{ sigma });
}

1;

package Statistics::Distribution::Generator::uniform;

use 5.010;
use base qw( Statistics::Distribution::Generator );
use overload (
  '0+' => '_render',
  '""' => '_render',
  '|' => '_add_alternative',
  'x' => '_add_dimension',
  fallback => 1,
);

sub _render {
  my $self = shift;
  return ($self->{ max } - $self->{ min }) * rand() + $self->{ min };
}

1;

package Statistics::Distribution::Generator::logistic;

use 5.010;
use base qw( Statistics::Distribution::Generator );
use overload (
  '0+' => '_render',
  '""' => '_render',
  '|' => '_add_alternative',
  'x' => '_add_dimension',
  fallback => 1,
);

sub _render {
  my $self = shift;
  return -log((1 / _rand_nonzero()) - 1);
}

1;

package Statistics::Distribution::Generator::supplied;

use 5.010;
use base qw( Statistics::Distribution::Generator );
use overload (
  '0+' => '_render',
  '""' => '_render',
  '|' => '_add_alternative',
  'x' => '_add_dimension',
  fallback => 1,
);

sub _render {
  my $self = shift;
  return $self->{ code }->();
}

1;

package Statistics::Distribution::Generator::gamma;

use 5.010;
use base qw( Statistics::Distribution::Generator );
use overload (
  '0+' => '_render',
  '""' => '_render',
  '|' => '_add_alternative',
  'x' => '_add_dimension',
  fallback => 1,
);

sub _gamma_int {
  my $order = shift;
  if ($order < 12){
    my $prod = 1;
    for (my $i=0; $i<$order; $i++){
      $prod *= _rand_nonzero();
    }
    return -log($prod);
  }
  else {
    return _gamma_large_int($order);
  }
}

sub _tan { sin($_[0]) / cos($_[0]); }

sub _gamma_large_int {
  my $order = shift;
  my $sqrt = sqrt(2 * $order - 1);
  my ($x,$y,$v);
  do {
    do {
      $y = _tan($pi * rand);
      $x = $sqrt * $y + $order - 1;
    } while ($x <= 0);
    $v = rand;
  } while ($v > (1 + $y * $y) * exp(($order - 1) * log($x / ($order - 1)) - $sqrt * $y));
  return $x;
}

sub _gamma_frac {
  my $order = shift;
  my $p = $e / ($order + $e);
  my ($q, $x, $u, $v);
  do {
    $u = rand;
    $v = _rand_nonzero();
    if ($u < $p){
      $x = exp((1 / $order) * log($v));
      $q = exp(-$x);
    }
    else {
      $x = 1 - log($v);
      $q = exp(($order - 1) * log($x));
    }
  } while (rand >= $q);
  return $x;
}

sub _render {
  my $self = shift;
  my $rv;
  if ($self->{ order } == $self->{ norder }) {
    $rv = $self->{ scale } * _gamma_int($self->{ norder });
  }
  elsif ($self->{ norder } == 0) {
    $rv = $self->{ scale } * _gamma_frac($self->{ order });
  }
  else {
    $rv = $self->{ scale } * (_gamma_int($self->{ norder }) + _gamma_frac($self->{ norder } - $self->{ order }));
  }
  return $rv;
}

1;

package Statistics::Distribution::Generator::exponential;

use 5.010;
use base qw( Statistics::Distribution::Generator );
use overload (
  '0+' => '_render',
  '""' => '_render',
  '|' => '_add_alternative',
  'x' => '_add_dimension',
  fallback => 1,
);

sub _render {
  my $self = shift;
  my $rv = -log(rand) / $self->{ lambda };
}

1;

package Statistics::Distribution::Generator::dice;

use 5.010;
use base qw( Statistics::Distribution::Generator );
use overload (
  '0+' => '_render',
  '""' => '_render',
  '|' => '_add_alternative',
  'x' => '_add_dimension',
  fallback => 1,
);

sub _render {
  my $self = shift;
  my $rv;
  for (1 .. $self->{ numdice }) {
    $rv += int(1 + rand($self->{ numsides }));
  }
  return $rv;
}

1;

__END__