|
The number sequence methods for generating random
numbers discussed in Chapter
4: Tech generally create uniform, or flat, distributions over
some range. However, what if you need some other type of non-uniform
distribution such as a Gaussian, Poisson, or even a custom distribution?
You can use one of two basic methods to generate a
non-uniform random distribution.
- Transformation Method - put the random value from the
uniform random value generator through a function to generate
a new value that follows a non-uniform distribution.
- Rejection Method - Best explained graphically (see below),
in this approach you need to generate two values y
and z
from unform random number generators. The first value y
is accepted as the output value if z
is less than p(y),
a standard curve that follows the desired shape. If z
is greater than p(y),
then y
is rejected and another pair of random values must be generated
and the test repeated.
Transformation
Method
The transformation method generate values of x
with a uniform random number generator and then transforms the x
values with a function y=g(x)
that results in a distribution following the desired probability
function f(y).
This method requires that the integral of f(y)exists
and that the integral be invertable. (The following derivation comes
from Press)
We define our probability function as p(y).
We want it to follow the shape of a particular function f(y)
such as a Gaussian or Poisson. That is,
p(y)
= f(y)
However, our uniform random number generator outputs
x with
p(x)
= 1
for x
between certain limits, such as 0 to 1, and is 0 otherwise. So we
need a function g(x)
y
= g(x)
that results in a distribution of y values that follows
the desired f(y).
Regardless of the transformation, the probabilities
will go as
|p(y)dy|
= |p(x)dx|
or
p(y)
= p(x) * dx/dy
Thus for the uniform random generator we obtain
p(y)
= dx/dy
= f(y)
An indefinite integral of f(y)*dy
gives
x = F(y)
Inverting this leads to
y(x)
= F-1 (x) = g(x)
So if our desired function is integrable and this
integral can be inverted, we have our analytical g(x)
transformation function. For example, for a sample of a radioactive
material, the number of decays would follow an exponential with
the half-life tau
p(y)
= exp(-t/tau) = f(y)
then
dx/dy
= exp(-y/tau)
This leads to
x
= tau(1 -
exp(-y/tau))
and then
y(x)
= -ln( 1 - x/tau)
for 0
<= x < tau. Since the (1
- x/tau) ranges randomly between 0 and 1.0, you can also
just use -ln(r)
for r generated
between 0 and 1.0.
The following applet shows a distribution generated with y(x)
= -ln(1-x/tau).
RanDistTransformApplet.java
- generate
a distribution in a histogram that follows a radioactive
decay distribution using the transform method discussed
above.
+ Previous
classes:
Chapter
6:Tech: Histogram.java,
HistPanel.java
Chapter
6:Tech: PlotPanel.java,
PlotFormat.java
|
|
import javax.swing.*;
import java.awt.*;
import java.awt.event.*;
/**
* This program will run as an applet
inside
* an application frame.
*
* The applet uses the HistPanel to display contents
of
* an instance of Histogram. HistFormat
used by HistPanel to
* format the scale values.
*
* Includes "Go" button to add random
values from a Gaussian
* distribution to the histogram.
The number of values taken from
* entry in a JTextField. "Clear" button
clears the histogram.
* In standalone mode, the Exit button
closes the program.
*
**/
public class RanDistTransformApplet extends JApplet
implements ActionListener
{
// Use the HistPanel JPanel subclass here
HistPanel fOutputPanel;
Histogram fHistogram ;
int fNumDataPoints = 1000;
// A text field for input strings
JTextField fTextField;
// Flag for whether the applet is in a browser
// or running via the main () below.
boolean fInBrowser = true;
//Buttons
JButton fGoButton;
JButton fClearButton;
JButton fExitButton;
/**
* Create a User Interface with histograms and
buttons to
* control the program. A textfield holds number
of entries
* to be generated for the histogram.
*/
public void init () {
Container content_pane = getContentPane
();
JPanel panel = new JPanel (new BorderLayout
());
// Create a histogram with Gaussian
distribution.
makeHist ();
// JPanel subclass here.
fOutputPanel = new HistPanel (fHistogram);
panel.add (fOutputPanel,"Center");
// Use a textfield for an input
parameter.
fTextField =
new JTextField (Integer.toString
(fNumDataPoints), 10);
// If return hit after entering
text, the
// actionPerformed will be invoked.
fTextField.addActionListener (this);
fGoButton = new JButton ("Go");
fGoButton.addActionListener (this);
fClearButton = new JButton ("Clear");
fClearButton.addActionListener (this);
fExitButton = new JButton ("Exit");
fExitButton.addActionListener (this);
JPanel control_panel = new JPanel
();
control_panel.add (fTextField);
control_panel.add (fGoButton);
control_panel.add (fClearButton);
control_panel.add (fExitButton);
if (fInBrowser) fExitButton.setEnabled
(false);
panel.add (control_panel,"South");
// Add text area with scrolling
to the content_pane.
content_pane.add (panel);
} // init
/** Respond to the buttons. **/
public void actionPerformed (ActionEvent e)
{
Object source = e.getSource ();
if (source == fGoButton || source
== fTextField) {
String strNumDataPoints
= fTextField.getText ();
try {
fNumDataPoints
= Integer.parseInt (strNumDataPoints);
}
catch (NumberFormatException
ex) {
//
Could open an error dialog here but just
//
display a message on the browser status line.
showStatus
("Bad input value");
return;
}
makeHist
();
repaint
();
} else if (source == fClearButton)
{
fHistogram.clear
();
repaint
();
} else if (!fInBrowser)
System.exit
(0);
} // actionPerformed
/** Create a histogram and fill it with data
using
* transformation algorithm to simulate
radioactive
* decay lifetimes.
**/
void makeHist () {
// Create an instance of the Random
class for
// producing our random values.
java.util.Random r = new java.util.Random
();
// Create an instance of our basic
histogram class.
// Make it wide enough enough to
include most of the
// gaussian values.
if (fHistogram == null)
fHistogram
= new Histogram ("Radioactive Lifetimes",
"Arbitrary
Units",
20,0.0,5.0);
// Use the transformation method
to generate a
// radioactive decay distribution
for (int i=0; i < fNumDataPoints;
i++) {
// y (x)
= -ln (1-x/tau);
// Generate
ran FP between 0 & 1
double x
= r.nextDouble ();
double val
= - Math.log (1-x);
fHistogram.add
(val);
}
} // makeHist
/** Use main() to run in standalone mode. **/
public static void main (String[] args) {
//
int frame_width=450;
int frame_height=300;
// Create an instance of this applet
an add to a frame.
RanDistTransformApplet applet =
new RanDistTransformApplet ();
applet.fInBrowser = false;
applet.init ();
// Following anonymous class used
to close window & exit program
JFrame f = new JFrame ("Demo");
f.setDefaultCloseOperation (JFrame.EXIT_ON_CLOSE);
// Add applet to the frame
f.getContentPane ().add ( applet);
f.setSize (new Dimension (frame_width,frame_height));
f.setVisible (true);
} // main
} // class RanDistTransformApplet
|
Note from the figure that instead of an analytical inversion, we
can use a graphical approach. Generate a uniform random value x
value along the vertical axis and then find the y that corresponds
to the integral up to that x value.
Gaussian Distributions
Though the java.util.Random
class provides the method nextGaussian()
to generate a random Gaussian distribution, it is still useful to
know how to write your own version in case you want to modify it
in some way.
The polar transformation [See Ref. Press and Tang] is the technique
used in java.util.Random
for generating Gaussians from uniform distributions. The Gaussian
function goes as
g(v)
= 1/sqrt(2*pi*sigma) * e[-v2/2sigma2]
A product of a uniform distribution p(x) over [0,1] with an exponential
distribution r(z)= e-z for z over the range [0, infinity]
can relate to a product of two Gaussians as follows:
p(x) dx r(z)dz = g(y1) dy1 g(y2)
dy2
This becomes (using sigma = 1.0):
2pi
* e-zdxdz = exp[-(y1)2 + (y2)2)/2]
dy1 dy2
Relating the two exponent factors, we see that this is equivalent
to a transformation from a rectangular coordinate system in
y1, y2
to a polar system where
r
= sqrt(2z)
and
y1
= sqrt(2z)cos(2pi*x)
y2 = sqrt(2z)sin(2pi*x)
The random z
exponential can be taken from an exponential distribution as shown
above:
y1
= sqrt(-2ln r)cos(2pi*x)
y2 =
sqrt(-2ln r)sin(2pi*x)
where now both r and x range from 0 to 1. So both y1
and y2
provide Gaussian distributed values.
A variation on this [Ref. Press] is to choose a random point (v1,v2)
within a unit circle. We use R = v12+v22 and then the cosine and
sine terms become v1/sqrt(R) and v2/sqrt(R), respectively. So
y1
= v1 * sqrt(-2(ln
R)/R)
y2 =
v2 * sqrt(-2(ln R)/R)
The applet here uses this technique:
RanDistGaussianApplet.java
-
generate a Gaussian distribution of numbers using a uniform
random number generator and the polar transform method
discussed above.
+ Previous classes:
Chapter
6:Tech: Histogram.java,
HistPanel.java
Chapter
6:Tech: PlotPanel.java,
PlotFormat.java
|
|
import javax.swing.*;
import java.awt.*;
import java.awt.event.*;
/**
* This program will run as an applet
inside
* an application frame.
*
* The applet uses the HistPanel to display contents
of
* an instance of Histogram. HistFormat
used by HistPanel to
* format the scale values.
*
* Includes "Go" button to add random
values from a Gaussian
* distribution to the histogram.
The number of values taken from
* entry in a JTextField. "Clear" button
clears the histogram.
* In standalone mode, the Exit button
closes the program.
*
**/
public class RanDistGaussianApplet extends JApplet
implements ActionListener
{
// Use the HistPanel JPanel subclass here
HistPanel fOutputPanel;
Histogram fHistogram;
int fNumDataPoints = 5000;
// A text field for input strings
JTextField fTextField;
// Flag for whether the applet is in a browser
// or running via the main () below.
boolean fInBrowser = true;
//Buttons
JButton fGoButton;
JButton fClearButton;
JButton fExitButton;
/**
* Create a User Interface with histograms and
buttons to
* control the program. A textfield holds number
of entries
* to be generated for the histogram.
*/
public void init () {
Container content_pane = getContentPane
();
JPanel panel = new JPanel (new BorderLayout
());
// Create a histogram with Gaussian
distribution.
makeHist ();
// JPanel subclass here.
fOutputPanel = new HistPanel (fHistogram);
panel.add (fOutputPanel,"Center");
// Use a textfield for an input
parameter.
fTextField =
new JTextField (Integer.toString
(fNumDataPoints), 10);
// If return hit after entering
text, the
// actionPerformed will be invoked.
fTextField.addActionListener (this);
fGoButton = new JButton ("Go");
fGoButton.addActionListener (this);
fClearButton = new JButton ("Clear");
fClearButton.addActionListener (this);
fExitButton = new JButton ("Exit");
fExitButton.addActionListener (this);
JPanel controlPanel = new JPanel
();
controlPanel.add (fTextField);
controlPanel.add (fGoButton);
controlPanel.add (fClearButton);
controlPanel.add (fExitButton);
if (fInBrowser) fExitButton.setEnabled
(false);
panel.add (controlPanel,"South");
// Add text area with scrolling
to the content_pane.
content_pane.add (panel);
}
/** Respond to the buttons. **/
public void actionPerformed (ActionEvent e)
{
Object source = e.getSource ();
if (source == fGoButton || source
== fTextField) {
String strNumDataPoints
= fTextField.getText ();
try {
fNumDataPoints
= Integer.parseInt (strNumDataPoints);
}
catch (NumberFormatException
ex) {
//
Could open an error dialog here but just
//
display a message on the browser status line.
showStatus
("Bad input value");
return;
}
makeHist
();
repaint
();
}
else if (source == fClearButton)
{
fHistogram.clear
();
repaint
();
} else if (!fInBrowser)
System.exit
(0);
} // actionPerformed
/** Creat the histogram and use Guassian - Polar
Transform
* algorithm to fill it with random
data.
**/
void makeHist () {
// Create an instance of the Random
class for
// producing our random values.
java.util.Random r = new java.util.Random
();
// Create an instance of our basic
histogram class.
// Make it wide enough enough to
include most of the
// gaussian values.
if (fHistogram == null)
fHistogram
= new Histogram ("Gaussian - Polar Transform",
"Arbitrary
Units",
25,-4.0,4.0);
// Use the transformation method
to generate a
// radioactive decay distribution
for (int i=0; i < fNumDataPoints;
i++) {
// Generate
random vals -1.0 to 1.0
double v1
= 2.0 * r.nextDouble () - 1.0;
double v2
= 2.0 * r.nextDouble () - 1.0;
// Restrict
them to inside a unit circle
double R
= v1*v1 + v2*v2;
if (R <
1.0) {
double
val = v1 * Math.sqrt (-2.0* (Math.log (R))/R);
fHistogram.add
(val);
}
}
} // makeHist
/** Use main() to run in standalone mode. **/
public static void main (String[] args) {
//
int frame_width=450;
int frame_height=300;
// Create an instance of this applet
an add to a frame.
RanDistGaussianApplet applet = new
RanDistGaussianApplet ();
applet.fInBrowser = false;
applet.init ();
// Following anonymous class used
to close window & exit program
JFrame f = new JFrame ("Demo");
f.setDefaultCloseOperation (JFrame.EXIT_ON_CLOSE);
// Add applet to the frame
f.getContentPane ().add ( applet);
f.setSize (new Dimension (frame_width,
frame_height));
f.setVisible (true);
} // main
} // class RanDistGaussianAppelt
|
Resources & Web References
Latest update: Nov. 4, 2004
|
