File Coverage

blib/lib/Statistics/GammaDistribution.pm

Criterion	Covered	Total	%
statement	71	75	94.6
branch	10	14	71.4
condition	3	9	33.3
subroutine	13	14	92.8
pod	4	6	66.6
total	101	118	85.5

package Statistics::GammaDistribution;
use strict;
use warnings;
use vars qw ( $VERSION );
$VERSION = 0.02;

sub PI(){ 3.14159265358979323846264338328; }
sub E() { 2.71828182845904523536028747135; }

sub new
{
    my ($caller,%args) = @_;
    my $class = ref($caller) || $caller;
    my $self = bless {}, $class;
    $self->_init(%args);
    return $self;
}

sub _init
{
    my ($self,%args) = @_;
    return $self;
}

sub get_order
{
    my ($self) = @_;
    return $self->{_a};
}

sub set_order
{
    my ($self,$order) = @_;
    $self->{_a} = $order;
    $self->{_na} = int($order);
    return $self->{_a};
}

sub rand
{
    my ($self,$scale) = shift;
    $scale = 1 unless defined $scale;
    my $order   = $self->{_a};
    my $n_order = $self->{_na};

    die('Statistics::GammaDistribution::rand() - order of distribution must be set greater than zero using set_order()')
	unless ((defined $order) && ($order>0));

    if ($order == $n_order){
	return $scale * _gamma_int($n_order);
    } elsif ($n_order == 0) {
	return $scale * _gamma_frac($order);
    } else {
	return $scale * (_gamma_int($n_order) + _gamma_frac($order-$n_order));
    }
}

# random number between zero and one,
# NOT including zero
sub _rand_nonzero
{
    my $rand;
    while(!($rand = CORE::rand()))
    {
	# loop
    }
    return $rand;
}

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 * CORE::rand());
	    $x = $sqrt * $y + $order - 1;
	} while ($x <= 0);
	$v = CORE::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 = CORE::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 (CORE::rand() >= $q);
    return $x;
}

sub dirichlet_dist
{
    my ($self,@alpha) = @_;
    my @theta;
    my $norm = 0.0;

    for(my $i=0; $i<=$#alpha; $i++){
	my $order = $alpha[$i];
	die('Statistics::GammaDistribution::dirichlet_dist() - every parameter must be greater than zero')
	    unless ((defined $order) && ($order>0));
	$self->set_order($order);
	$theta[$i] = $self->rand(1);
	$norm += $theta[$i];
    }

    for(my $i=0; $i<=$#theta; $i++){
	$theta[$i] /= $norm;
    }

    return @theta;
}

1;

__END__