Bio/Tools/RandomDistFunctions.pm

   1 # $Id$
   2 #
   3 # BioPerl module for Bio::Tools::RandomDistFunctions
   4 #
   5 # Cared for by Jason Stajich <jason-at-bioperl.org>
   6 #
   7 # Copyright Jason Stajich
   8 #
   9 # You may distribute this module under the same terms as perl itself
  10
  11 # POD documentation - main docs before the code
  12
  13 =head1 NAME
  14
  15 Bio::Tools::RandomDistFunctions - A set of routines useful for
  16 generating random data in different distributions
  17
  18 =head1 SYNOPSIS
  19
  20   use Bio::Tools::RandomDistFunctions;
  21   my $dist = Bio::Tools::RandomDistFunctions->new();
  22   for my $v ( 1..1000 ) {
  23     my $birth_dist = $dist->rand_birth_distribution($lambda);
  24     # ... do something with the variable
  25   }
  26
  27 =head1 DESCRIPTION
  28
  29 Most of the code is based on the C implementation of these routines in
  30 Mike Sanderson's r8s's package.  See http://ginger.ucdavis.edu/ for
  31 information on his software.
  32
  33 This code tries to be fast and use available faster BigInt and GMP
  34 library methods when those modules are available.
  35
  36 =head1 FEEDBACK
  37
  38 =head2 Mailing Lists
  39
  40 User feedback is an integral part of the evolution of this and other
  41 Bioperl modules. Send your comments and suggestions preferably to
  42 the Bioperl mailing list.  Your participation is much appreciated.
  43
  44   bioperl-l@bioperl.org                  - General discussion
  45   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  46
  47 =head2 Reporting Bugs
  48
  49 Report bugs to the Bioperl bug tracking system to help us keep track
  50 of the bugs and their resolution. Bug reports can be submitted via the
  51 web:
  52
  53   http://bugzilla.open-bio.org/
  54
  55 =head1 AUTHOR - Jason Stajich
  56
  57 Email jason-at-bioperl.org
  58
  59 =head1 CONTRIBUTORS
  60
  61 Thanks to Mike Sanderson for assistance in the getting this
  62 implementation together.
  63
  64 =head1 APPENDIX
  65
  66 The rest of the documentation details each of the object methods.
  67 Internal methods are usually preceded with a _
  68
  69 =cut
  70
  71
  72 # Let the code begin...
  73
  74
  75 package Bio::Tools::RandomDistFunctions;
  76 require Exporter;
  77 use vars qw(%LOADED @EXPORT_OK); use strict;
  78
  79 #use Math::BigFloat lib => 'GMP,Bit::Vector';
  80 #use Math::BigInt lib  => 'GMP,Bit::Vector';
  81 use POSIX;
  82
  83 use base qw(Bio::Root::Root);
  84
  85 =head2 birth_distribution
  86
  87  Title   : rand_birth_distribution
  88  Usage   : my $randvar = $dist->
  89             rand_birth_distribution($lambda);
  90  Function: Returns a random number from a birth process waiting
  91            time with a fixed interval
  92            1.0.  Times are measured from 0=present,1=root;
  93  Returns : floating point number
  94  Args    : $lambda ( > 0 )
  95  References : This is based on code by Mike Sanders in r8s.
  96               Ross, Stochastic Processes, p. 145 for the density
  97
  98 =cut
  99
 100 sub rand_birth_distribution{
 101    my ($self,$lambda) = @_;
 102    if( ! ref($self) &&
 103        $self !~ /RandomDistFunctions/ ) {
 104        $lambda = $self;
 105    }
 106    unless( $lambda ) {
 107        $self->throw("Cannot call birth_distribution without a valid lambda value (>0)");
 108    }
 109    return 1 - (log(rand(1) * (exp($lambda) - 1)+1)/ $lambda);
 110 }
 111
 112
 113 =head2 rand_geometric_distribution
 114
 115  Title   : rand_geometric_distribution
 116  Usage   : my $randvar = $dist->rand_geometric_distribution($param);
 117  Function: Returns a random geometric variate distributed with
 118            paramater $param, according to
 119            c.d.f. 1 - ( 1- param) ^ n
 120  Returns : integer
 121  Args    : $param ( 0 > $param < 1 )
 122
 123
 124 =cut
 125
 126 sub rand_geometric_distribution{
 127    my ($self,$param) = @_;
 128    if( ! ref($self) &&
 129        $self !~ /RandomDistFunctions/ ) {
 130        $param = $self;
 131    }
 132    unless( $param ) {
 133        $self->throw("Cannot call rand_geometric_distribution without a valid param value (>0)");
 134    }
 135
 136    my $den;
 137    if( $param < 1e-8) {
 138        $den = (-1 * $param) - ( $param * $param ) / 2;
 139    } else {
 140        $den = log(1 - $param);
 141    }
 142    my $z = log(1 - rand(1)) / $den;
 143    return POSIX::floor($z) + 1;
 144    # MSanderson comments from r8s code
 145    # Is this the right truncation of the real-valued expression above?
 146    # YES
 147    # Checked by reference to the expected mean of the distribution in
 148    # 100,000 replicates
 149    # EX = 1/param Var = (1-param)/param^2 See Olkin, Gleser, and
 150    # Derman, p. 193ff. Probability Models and Applications, 1980.
 151 }
 152
 153 =head2 rand_exponentional_distribution
 154
 155  Title   : rand_exponentional_distribution
 156  Usage   : my $var = $dist->rand_exponentional_distribution($param);
 157  Function: Returns a random exponential variate distributed with parameter
 158            $param, according to c.d.f 1 - e^(-param * x)
 159  Returns : floating point number
 160  Args    : $param ( > 0 )
 161
 162
 163 =cut
 164
 165 sub rand_exponentional_distribution {
 166    my ($self,$param) = @_;
 167       if( ! ref($self) &&
 168        $self !~ /RandomDistFunctions/ ) {
 169        $param = $self;
 170    }
 171    unless( $param ) {
 172        $self->throw("Cannot call rand_exponentional_distribution without a valid param value (>0)");
 173    }
 174    return log( 1- rand(1)) / $param;
 175 }
 176
 177 =head2 rand_normal_distribution
 178
 179  Title   : rand_normal_distribution
 180  Usage   : my $var = $dist->rand_normal_distribution()
 181  Function: Returns a random normal (gaussian) variate distributed
 182  Returns : floating point number
 183  Args    : none
 184
 185
 186 =cut
 187
 188 sub rand_normal_distribution{
 189     my $gset;
 190     my ($rsq,$v1,$v2) = ( 0,0,0);
 191     do {
 192         $v1 = 2 * rand(1) - 1;
 193         $v2 = 2 * rand(1) - 1;
 194         $rsq= $v1**2 + $v2 ** 2;
 195     } while( $rsq >= 1 || $rsq == 0);
 196     my $fac = sqrt(-2.0 * log($rsq) / $rsq );
 197     $gset = $v1 * $fac;
 198     return $v2 * $fac;
 199 }
 200
 201 1;