Statistics::GammaDistribution - represents a gamma distribution

Statistics-GammaDistribution documentation Contained in the Statistics-GammaDistribution distribution.

Index

Code Index:

NAME

Top

Statistics::GammaDistribution - represents a gamma distribution

SYNOPSIS

Top

  use Statistics::GammaDistribution;
  my $g = Statistics::GammaDistribution->new();
  $g->set_order(8.5);
  print $g->rand(1.0);

  my @alpha = (0.5,4.5,20.5,6.5,1.5,0.5);
  my @theta = $g->dirichlet_dist(@alpha);

METHODS

Top

$gamma = Statistics::GammaDistribution->new();

No parameters necessary.

$variate = $gamma->rand( SCALE );

This function returns a random variate from the gamma distribution. The distribution function is,

  p(x) dx = {1 \over \Gamma(a) b^a} x^{a-1} e^{-x/b} dx
  for x > 0.

Where a is the order and b is the scale. Unless supplied as a parameter, SCALE is assumed to be 1.0 if not supplied.

$gamma->get/set_order( ORDER );

Gets/sets the order of the distribution. Order must be greater than zero.

@theta = $gamma->dirichlet_dist( ALPHA );

Takes a K-sized array of real numbers (all greater than zero), and returns a K-sized array containing random variates from a Dirichlet distribution. The distribution function is

  p(\theta_1, ..., \theta_K) d\theta_1 ... d\theta_K = 
    (1/Z) \prod_{i=1}^K \theta_i^{\alpha_i - 1} \delta(1 -\sum_{i=1}^K \theta_i) d\theta_1 ... d\theta_K

    for theta_i >= 0 and alpha_i >= 0. The normalization factor Z is
  Z = {\prod_{i=1}^K \Gamma(\alpha_i)} / {\Gamma( \sum_{i=1}^K \alpha_i)}

The random variates are generated by sampling K values from gamma distributions with parameters order=alpha_i, scale=1, and renormalizing. See A.M. Law, W.D. Kelton, Simulation Modeling and Analysis (1991).

AUTHOR

Top

Nigel Wetters Gourlay <nwetters@cpan.org>

COPYRIGHT

Top

Copyright (c) 2003-06, Nigel Wetters Gourlay. NO WARRANTY.

This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

Statistics-GammaDistribution documentation Contained in the Statistics-GammaDistribution distribution. 
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__