nr02doc.htm

随机数发生器C++写的

  </tr>
  <tr>
    <td width="50%" valign="top">Microsoft 5.0, 6.0, 7.0</td>
    <td width="50%" valign="top">OK</td>
  </tr>
  <tr>
    <td width="50%" valign="top">Watcom 10A</td>
    <td width="50%" valign="top">OK</td>
  </tr>
  <tr>
    <td width="50%" valign="top">Sun CC 6</td>
    <td width="50%" valign="top">OK</td>
  </tr>
  <tr>
    <td width="50%" valign="top">Gnu G++ 2.96 (Linux), 2.95 (Sun)</td>
    <td width="50%" valign="top">OK</td>
  </tr>
  <tr>
    <td width="50%" valign="top">Intel 5 for Windows or Linux</td>
    <td width="50%" valign="top">OK</td>
  </tr>
</table>

<p>I have included make files for Borland 5.5, CC and Gnu G++. See <a href="#files">files</a>
section. You can use my <a href="genmake.htm">genmake</a> program to generate 
make files for other compilers.</p>

<h3><a name="testing"></a>Testing</h3>

<p>The files tryrand.cpp, tryrand1.cpp, tryrand2.cpp, tryrand3.cpp, 
tryrand4.cpp, tryrand5.cpp, hist.cpp
run the generators in the library and print histograms of the resulting distributions.
Sample means and variances are also calculated and can be compared with the population
values. The  results from Borland 5.5 are in tryrand.txt. Other compilers may 
give different but still correct results. In particular, Microsoft C++ compilers 
give different results when global optimisation is set. This is most likely due 
to different round-off error but may also be due to different order of 
evaluation of expressions. The appearance of the histograms in the output should 
be similar to that in tryrand.txt and the statistical tests should still be 
passed. There are some notes in tryrand.txt for assessing the output.&nbsp;
&nbsp; </p>

<table BORDER="1" WIDTH="100%">
  <tr>
    <td WIDTH="100%">If you are compiling on a PC with a 16 bit compiler you will need to set <i>n_large</i>
    to be 8000 rather than 1000000 and <i>n</i> to 8000 rather than 200000 in tryrand.cpp.
    This will change the output from the test program from that given in tryrand.txt although
    the general appearance should be the same.</td>
  </tr>
</table>

<p>The test program tryrand.cpp includes a simple test for memory leaks. This is valid for
only some compilers. It seems to work for Borland C++ in console mode but not for Gnu G++
or Microsoft C++, where it almost always (presumably incorrectly) suggests an error. </p>

<p>&nbsp; </p>

<h2><a NAME="description"></a>Descriptions of the classes to be
accessed by the user:</h2>

<h3><a NAME="random"></a>Random:</h3>

<p>This is the basic uniform random number generator, used to drive all the others. The
Lewis-Goodman-Miller algorithm is used with Marsaglia mixing. While not perfect, and now
superseded, the LGM generator has given me acceptable results in a wide variety of
simulations. See <i>Numerical Recipes in C</i> by Press, Flannery, Teukolsky, Vetterling
published by the Cambridge University Press for details. The LGM generator does pass the
Marsaglia <i><a HREF="http://stat.fsu.edu/~geo/diehard.html">diehard</a></i> tests when
you include the mixing. (It doesn't pass without the mixing). Nevertheless it should be
upgraded. It is very dubious if you are going to call it more than 100 million 
times in a simulation. Ideally the basic generator should be recoded in assembly language to give the
maximum speed to all the generators in this package. You can access the numbers directly
using <tt>Next()</tt> but I suggest you use class <a HREF="#uniform">Uniform</a> for
uniform random numbers and reserve Random for setting the starting seed and as the base
class for the random number generators.</p>

<p>Remember that you need to call <tt>Random::Set(double seed)</tt> at the beginning of 
your program. See the <a href="#overview">overview</a> and tryrand.cpp.</p>

<h3><a NAME="uniform"></a>Uniform:</h3>

<p>Return a uniform random number from the range (0, 1). The constructor has no
parameters. For example </p>

<pre>&nbsp;&nbsp; Uniform U;
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; U.Next() &lt;&lt; &quot;\n&quot;;</pre>

<p>prints a column of 100 numbers drawn from a uniform distribution. </p>

<h3><a NAME="constant"></a>Constant:</h3>

<p>This returns a constant. The constructor takes one <i>Real</i> parameter; the value of
the constant to be returned. So </p>

<pre>&nbsp;&nbsp; Constant C(5.5);
&nbsp;&nbsp; cout &lt;&lt; C.Next() &lt;&lt; &quot;\n&quot;;</pre>

<p>prints 5.5. </p>

<h3><a NAME="exponential"></a>Exponential:</h3>

<p>This generates random numbers with density <tt>exp(-x)</tt> for <tt>x&gt;=0</tt>. The
constructor takes no arguments. </p>

<pre>&nbsp;&nbsp; Exponential E;
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; E.Next() &lt;&lt; &quot;\n&quot;;</pre>

<h3><a NAME="cauchy"></a>Cauchy:</h3>

<p>Generates random numbers from a standard Cauchy distribution. The constructor takes no
parameters. </p>

<pre>&nbsp;&nbsp; Cauchy C;
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; C.Next() &lt;&lt; &quot;\n&quot;;</pre>

<h3><a NAME="normal"></a>Normal:</h3>

<p>Generates standard normal random numbers. The constructor has no arguments. This class
has been augmented to ensure only one copy of the arrays generated by the constructor
exist at any given time. That is, if the constructor is called twice (before the
destructor is called) only one copy of the arrays is generated. </p>

<pre>&nbsp;&nbsp; Normal Z;
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; Z.Next() &lt;&lt; &quot;\n&quot;;</pre>

<h3><a NAME="chisq"></a>ChiSq:</h3>

<p>Non-Central chi-squared distribution. The method uses ChiSq1 to generate the
non-central part and Gamma2 or Exponential to generate the central part. The constructor
takes as arguments the number of degrees of freedom <tt>(&gt;=1)</tt> and the
non-centrality parameter (omit if zero). </p>

<pre>&nbsp;&nbsp; int df = 10; Real noncen = 2.0;
&nbsp;&nbsp; ChiSq CS(df, noncen);
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; CS.Next() &lt;&lt; &quot;\n&quot;;</pre>

<h3><a NAME="gamma"></a>Gamma:</h3>

<p>Gamma distribution. The constructor takes the shape parameter as argument. Uses Gamma1,
Gamma2 or Exponential. </p>

<pre>&nbsp;&nbsp; Real shape = 0.75;
&nbsp;&nbsp; Gamma G(shape);
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; G.Next() &lt;&lt; &quot;\n&quot;;</pre>

<h3><a NAME="pareto"></a>Pareto:</h3>

<p>Pareto distribution. The constructor takes the shape parameter as argument. I follow
the definition of Kotz and Johnson's <i>Continuous univariate distributions 1</i>, chapter
19, page 234, with <i>k</i> = 1. The generator uses a power transform of a uniform random
number. </p>

<pre>&nbsp;&nbsp; Real shape = 0.75;
&nbsp;&nbsp; Pareto P(shape);
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; P.Next() &lt;&lt; &quot;\n&quot;;</pre>

<h3><a NAME="poisson"></a>Poisson:</h3>

<p>Poisson distribution: uses Poisson1 or Poisson2. Constructor takes the mean as its
argument. </p>

<pre>&nbsp;&nbsp; Real mean = 5.0;
&nbsp;&nbsp; Poisson P(mean);
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; (int)P.Next() &lt;&lt; &quot;\n&quot;;</pre>

<h3><a NAME="binomial"></a>Binomial:</h3>

<p>Binomial distribution: uses Binomial1 or Binomial2. Constructor takes <i>n</i> and <i>p</i>
as its arguments. </p>

<pre>&nbsp;&nbsp; int n = 50; Real p = 0.25;
&nbsp;&nbsp; Binomial B(n, p);
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; (int)B.Next() &lt;&lt; &quot;\n&quot;;</pre>

<h3><a NAME="negativebinomial"></a>NegativeBinomial:</h3>

<p>Negative binomial distribution. Constructor takes <i>N</i> and <i>P</i> as its
arguments. I use the notation of Kotz and Johnson's <i>Discrete distributions</i>. Some
people use <i>p</i> = 1/(<i>P</i>+1) in place of the second parameter. </p>

<pre>&nbsp;&nbsp; Real N = 12.5; Real P = 3.0;
&nbsp;&nbsp; NegativeBinomial NB(N, P);
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; (int)NB.Next() &lt;&lt; &quot;\n&quot;;</pre>

<h3><a NAME="posgenx"></a>PosGenX:</h3>

<p>This uses an arbitrary density satisfying the previous conditions to generate random
numbers from that density. Suppose <tt>Real pdf(Real)</tt> is the density. Then use <tt>pdf</tt>
as the argument of the constructor. For example </p>

<pre>&nbsp;&nbsp; PosGenX P(pdf);
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; P.Next() &lt;&lt; &quot;\n&quot;;</pre>

<table BORDER="1" WIDTH="100%">
  <tr>
    <td WIDTH="100%">Note that the probability density <i>pdf</i> must drop to exactly 0 for
    the argument large enough. For example, include a statement in the program for <i>pdf</i>
    that, if the value is less than 1.0E-15, then return 0.</td>
  </tr>
</table>

<h3><a NAME="symgenx"></a>SymGenX:</h3>

<p>This corresponds to PosGenX for symmetric distributions. <br>
&nbsp; </p>

<table BORDER="1" WIDTH="100%">
  <tr>
    <td WIDTH="100%">Note that the probability density <i>pdf</i> must drop to exactly 0 for
    the argument large enough. For example, include a statement in the program for <i>pdf</i>
    that, if the value is less than 1.0E-15, then return 0.</td>
  </tr>
</table>

<h3><a NAME="asymgenx"></a>AsymGenX:</h3>

<p>Corresponds to PosGenX. The arguments of the constructor are the name of the density
function and the location of the mode. </p>

<pre>&nbsp;&nbsp; Real pdf(Real);
&nbsp;&nbsp; Real mode;
&nbsp;&nbsp; .....
&nbsp;&nbsp; AsymGenX X(pdf, mode);
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; X.Next() &lt;&lt; &quot;\n&quot;;</pre>

<table BORDER="1" WIDTH="100%">
  <tr>
    <td WIDTH="100%">Note that the probability density <i>pdf</i> must drop to exactly 0 for
    the argument large (large positive and large negative) enough. For example, include a
    statement in the program for <i>pdf</i> that, if the value is less than 1.0E-15, then
    return 0.</td>
  </tr>
</table>

<h3><a NAME="posgen"></a>PosGen:</h3>

<p>PosGen is not used directly. It is used as a base class for generating a random number
from an arbitrary probability density <tt>p(x)</tt>. <tt>p(x)</tt> must be non-zero only
for <tt>x&gt;=0</tt>, be monotonically decreasing for <tt>x&gt;=0</tt>, and be finite. For
example, <tt>p(x)</tt> could be <tt>exp(-x)</tt> for <tt>x&gt;=0</tt>. </p>

<p>The method is to cover the density in a set of rectangles of equal area as in the
diagram (indicated by <tt>---</tt>). </p>

<pre>&nbsp;&nbsp; <b>|
</b>&nbsp;&nbsp; <b>x
</b>&nbsp;&nbsp; <b>|xx</b>------
&nbsp;&nbsp; <b>|</b>&nbsp; <b>xx</b>&nbsp;&nbsp;&nbsp; |
&nbsp;&nbsp; <b>|</b>&nbsp;&nbsp;&nbsp; <b>xxx</b> |
&nbsp;&nbsp; <b>|</b>.......<b>xxx</b>---------
&nbsp;&nbsp; <b>|</b>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; | <b>xxxx</b>&nbsp;&nbsp;&nbsp;&nbsp; |
&nbsp;&nbsp; <b>|</b>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp; <b>xxxx</b> |
&nbsp;&nbsp; <b>|</b>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |.........<b>xxxxx</b>------------
&nbsp;&nbsp; <b>|</b>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp; <b>xxxxx</b>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |
&nbsp;&nbsp; <b>|</b>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <b>xxxxxx</b> |
&nbsp;&nbsp; <b>|</b>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |..............<b>xxxxxx</b>----------------------
&nbsp;&nbsp; <b>|</b>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp; <b>xxxxxxx</b>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |
&nbsp;&nbsp; <b>|</b>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <b>xxxxxxx</b>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |
&nbsp;&nbsp; <b>|</b>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; |&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <b>xxxxxxxx</b>|
&nbsp;&nbsp; <b>+===========================================================================</b></pre>

<p>The numbers are generated by generating a pair of numbers uniformly distributed over
these rectangles and then accepting the <i>X</i> coordinate as the next random number if
the pair corresponds to a point below the density function. The acceptance can be done in
two stages, the first being whether the number is below the dotted line. This means that
the density function need be checked only very occasionally and on the average only just
over 3 uniform random numbers are required for each of the random numbers produced by this
generator. </p>

<p>See PosGenX or Exponential for the method of deriving a class to generate random
numbers from a given distribution. <br>
&nbsp; </p>

<table BORDER="1" WIDTH="100%">
  <tr>
    <td WIDTH="100%">Note that the probability density <i>p(x)</i> must drop to exactly 0 for
    the argument, <i>x</i>, large enough. For example, include a statement in the program for <i>p(x)</i>
    that, if the value is less than 1.0E-15, then return 0.</td>
  </tr>
</table>

<h3><a NAME="symgen"></a>SymGen:</h3>