Re: Ask Dr. Random

[Paul E Johnson]

I would second Marcus's recommendation to look at the GSL. For a while,
I was writing up random number classes to donate to Swarm,but I could
not work my way through all those details (that Sven and Marcus manage
on a daily basis).  It is pretty easy to get it setup to draw you a
random number. The code snip below is from a project where I create a
Beta distribution and a Binomial distribution. I took the function from
the GSL project for Binomial and it was easy to craft it to an ObjC
method.

If you don't have the time to do something like this for Poisson,
contact me (address@hidden) and I'll do it for you if you promise to
put my picture on your web page.:)


Paul Johnson

- initializeRNGs
{
  if(randomSeed != -1)
    {
    [randomGenerator setStateFromSeed: randomSeed];
    }
  else
    {
    randomSeed= [randomGenerator getInitialSeed];
    
    //printf("Random number seed is builtin %d \n", randomSeed);
    }
  GammaDist1= [GammaDist create: self  setGenerator: randomGenerator
setAlpha: 2 setBeta: 1]; 
  return self;
}



-(double) getBetaVariateAlpha1: (double) alpha1 Alpha2: (double) alpha2
{
  double y1, y2;
  y1=[GammaDist1 getSampleWithAlpha: alpha1  withBeta: 1];
  y2=[GammaDist1 getSampleWithAlpha: alpha2  withBeta: 1];
  return y1/(y1+y2);    
}


//Paul Johnson July 13,2000.  This next method "getBinomialVariateN:
//P: is adapted from the C code of the GSL, which is available under
//the GPL with this statement:

/* randist/binomial.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
 *
 * 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.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* The binomial distribution has the form,
   prob(k) =  n!/(k!(n-k)!) *  p^k (1-p)^(n-k) for k = 0, 1, ..., n
   This is the algorithm from Knuth */


-(int) getBinomialVariateN: (int)  n P: (double) p
{
  unsigned int i, a, b, k = 0;
  while (n > 10)        /* This parameter is tunable */
    {
      double X;
      a = 1 + (n / 2);
      b = 1 + n - a;
      X = [self getBetaVariateAlpha1: a Alpha2: b];
      if (X >= p)
        {
          n = a - 1;
          p /= X;
        }
      else
        {
          k += a;
          n = b - 1;
          p = (p - X) / (1 - X);
        }
    }
  for (i = 0; i < n; i++)
    {
      double u = [uniformDblRand getDoubleWithMin: 0 withMax: 1];
      if (u < p)
        k++;
    }
  return k;
}


Miles Parker wrote:
> 
> Let me be the first to say that I think the Poisson distribution would be a
> perfect choice, even if it does sound a bit fishy!
> 
> Miles T. Parker


-- 
Paul E. Johnson                       email: address@hidden
Dept. of Political Science            http://lark.cc.ukans.edu/~pauljohn
University of Kansas                  Office: (785) 864-9086
Lawrence, Kansas 66045                FAX: (785) 864-5700


                  ==================================
   Swarm-Modelling is for discussion of Simulation and Modelling techniques
   esp. using Swarm.  For list administration needs (esp. [un]subscribing),
   please send a message to <address@hidden> with "help" in the
   body of the message.
                  ==================================