Login | Register   
LinkedIn
Google+
Twitter
RSS Feed
Download our iPhone app
TODAY'S HEADLINES  |   ARTICLE ARCHIVE  |   FORUMS  |   TIP BANK
Browse DevX
Sign up for e-mail newsletters from DevX


Tip of the Day
Language: Math processing
Expertise: Beginner
Mar 19, 1997

Including Random Numbers from Poisson in Application

Question:
How can I include random numbers drawn from a Poisson distribution in my applications? There is only a Gaussian distribution in java.util.Random.

Answer:
Well, this one is a little off the beaten track. I don't know how many people routinely need to generate Poisson distributed random numbers, so it isn't too surprising that there is no standard function for it. (The Poisson distribution describes probabilities of certain rare events. The French mathematician Poisson developed the mathematical theory when studying the frequency of soldiers being killed by horse kicks to their heads in the Napoleonic army.)

An excellent reference for this kind of problems is "Numerical Recipes in C", by Press, Teukolsky, Vetterling and Flannery, published by Cambridge University Press.

In the sample program, I implemented a class Poisson. To make Poisson distributed numbers, simply make an object of this class, giving the mean in the constructor:

        Poisson p = new Poisson(5);
        

Then keep calling

        p.next()
        

to get the random numbers.

Here is an applet that gives a visual idea of the Poisson distribution. Enter a mean (such as 5, 20, or 50) in the "Mean field", click on "Generate" and wait a short while. You will get a histogram of 1,000 trials, showing how many times the numbers 1, 2, 3, ..., 100 were generated by the Poisson process.

import java.applet.*;
import java.awt.*;

public class PoissonApplet extends Applet
{  public void init()
   {  graph = new HistoCanvas();
      setLayout(new BorderLayout());
      Panel p = new Panel();
      mean = new TextField("20", 10);
      p.add(new Label("Mean: "));
      p.add(mean);
      p.add(new Button("Generate"));
      add("North", p);
      add("Center", graph);
   }

   public void start()
   {  generate();
   }

   public boolean action(Event evt, Object arg)
   {  if (arg.equals("Generate"))
      {  generate();
      }
      else return super.action(evt, arg);
      return true;
   }

   public void generate()
   {  int[] v = new int[100];
      try
      {  float m = Float.valueOf(mean.getText().trim()).floatValue();
         Poisson p1 = new Poisson(m);
         for (int i = 0; i < 2000; i++)
         {  int p = p1.next();
            if (0 <= p && p < v.length) v[p]++;
         }
         graph.redraw(v, 0, 500);
      } catch (NumberFormatException e) {}
   }

   HistoCanvas graph;
   TextField mean;
}

class Poisson
{  public Poisson(float m)
   {  mean = m;
      if (mean < 12)
      {  g = (float)Math.exp(-mean);
      }
      else
      {  sq = (float)Math.sqrt(2 * mean);
         lm = (float)Math.log(mean);
         g = mean * lm - loggamma(mean + 1);
      }
   }

   public int next()
   {  if (mean < 12)
      {  int e = -1;
         float t = 1;
         do
         {  e++;
            t *= Math.random();
         } while (t > g);
         return e;
      }
      else
      {  float e;
         float t;
         do
         {  float y;
            do
            {  y = (float)Math.tan(Math.PI * Math.random());
               e = sq * y + mean;
            } while (e < 0);
           e = (float)Math.floor(e);
            t = 0.9F * (1.0F + y * y) * (float)Math.exp(e * lm - 
loggamma(e
 1) - g);
         } while (Math.random() > t);
         return (int)e;
      }
   }


   private static float loggamma(float xx)
   {  double x = xx;
      double y = x;
      double tmp = x + 5.5;
      tmp -= (x + 0.5) * Math.log(tmp);
      double ser = 1.000000000190015;
      for (int i = 0; i < 6; i++)
      {  y++;
         ser += coeff[i] / y;
      }
      return (float)(-tmp + Math.log(2.5066282746310005 * ser / x));
   }

   private float mean;
   private float g;
   private float sq;
   private float lm;

   private static double[] coeff;
   static
   {  coeff = new double[6];
      coeff[0] = 76.18009172947146;
      coeff[1] = -86.50532032941677;
      coeff[2] = 24.01409824083091;
      coeff[3] = -1.231739572450155;
      coeff[4] = 0.1208650973866179E-2;
      coeff[5] = -0.5395239384953E-5;
   }

   public static void main(String[] args)
   {  // test loggamma
      System.out.println(loggamma(1));
      System.out.println(loggamma(10));
      System.out.println(Math.log(2*3*4*5*6*7*8*9));
      // test Poisson
      Poisson p1 = new Poisson(5);
      float sum = 0;
      for (int i = 0; i < 10; i++) System.out.println(p1.next());
      for (int i = 0; i < 100; i++) sum += p1.next();
      System.out.println(sum / 100);
      Poisson p2 = new Poisson(15);
      sum = 0;
      for (int i = 0; i < 10; i++) System.out.println(p2.next());
      for (int i = 0; i < 100; i++) sum += p2.next();
      System.out.println(sum / 100);
   }
}

class HistoCanvas extends Canvas
{  public HistoCanvas()
   {  resize(200, 200);
   }

   public void redraw(int[] values, int ymin, int ymax)
   {  info = values;
      minValue = ymin;
      maxValue = ymax;
      repaint();
   }

   public void paint(Graphics g)
   {  if (info == null) return;

      int i;
      Dimension d = size();
      int barWidth = d.width / info.length;
      double scale = (double)d.height
         / (maxValue - minValue);

      for (i = 0; i < info.length; i++)
      {  int x1 = i * barWidth + 1;
         int y1;
         int v = info[i];
         int height;
         int yOrigin = (int)(maxValue * scale);

         if (v >= 0)
         {  y1 = (int)((maxValue - v) * scale);
            height = yOrigin - y1;
         }
         else
         {  y1 = yOrigin;
            height = (int)(-v * scale);
         }

         g.setColor(Color.red);
         g.fillRect(x1, y1, barWidth - 2, height);
         g.setColor(Color.black);
         g.drawRect(x1, y1, barWidth - 2, height);
      }
   }

   private int[] info = null;
   private int minValue;
   private int maxValue;
}
DevX Pro
 
Comment and Contribute

 

 

 

 

 


(Maximum characters: 1200). You have 1200 characters left.

 

 

Sitemap
Thanks for your registration, follow us on our social networks to keep up-to-date