Bio/Tools/RandomDistFunctions.pm

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