nr02doc.htm

随机数生成程序

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

<h3><a NAME="constant"></a><font COLOR="#FF0000">Constant:</font></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><font COLOR="#FF0000">Exponential:</font></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><font COLOR="#FF0000">Cauchy:</font></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><font COLOR="#FF0000">Normal:</font></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><font COLOR="#FF0000">ChiSq:</font></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><font COLOR="#FF0000">Gamma:</font></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><font COLOR="#FF0000">Pareto:</font></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><font COLOR="#FF0000">Poisson:</font></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><font COLOR="#FF0000">Binomial:</font></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><font COLOR="#FF0000">NegativeBinomial:</font></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><font COLOR="#FF0000">PosGenX:</font></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><font COLOR="#FF0000">SymGenX:</font></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><font COLOR="#FF0000">AsymGenX:</font></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><font COLOR="#FF0000">PosGen:</font></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><font COLOR="#FF0000">SymGen:</font></h3>

<p>SymGen is a modification of PosGen for unimodal distributions symmetric about the
origin, such as the standard normal. </p>

<h3><a NAME="asymgen"></a><font COLOR="#FF0000">AsymGen:</font></h3>

<p>A general random number generator for unimodal distributions following the method used
by PosGen. The constructor takes one argument: the location of the mode of the
distribution. </p>

<h3><a NAME="discretegen"></a><font COLOR="#FF0000">DiscreteGen:</font></h3>

<p>This is for generating random numbers taking just a finite number of values. There are
two alternative forms of the constructor: </p>

<pre>&nbsp;&nbsp; DiscreteGen D(n,prob);
&nbsp;&nbsp; DiscreteGen D(n,prob,val);</pre>

<p>where <tt>n</tt> is an integer giving the number of values, <tt>prob</tt> a Real array
of length <tt>n</tt> giving the probabilities and <tt>val</tt> a Real array of length <tt>n</tt>
giving the set of values that are generated. If <tt>val</tt> is omitted the values are <tt>0,1,...,n-1</tt>.
</p>

<p>The method requires two uniform random numbers for each number it produces. This method
is described by Kronmal and Peterson, <i>American Statistician</i>, 1979, Vol 33, No 4,
pp214-218. </p>

<h3><a NAME="sum"></a><font COLOR="#FF0000">SumRandom:</font></h3>

<p>This is for building a random number generator as a linear or multiplicative
combination of existing random number generators. Suppose <tt>RV1</tt>, <tt>RV2</tt>, <tt>RV3</tt>,
<tt>RV4</tt> are random number generators defined with constructors given above and <tt>r1</tt>,
<tt>r2</tt>, <tt>r0</tt> are Reals and <tt>i1</tt>, <tt>i3</tt> are integers. </p>

<p>Then the generator <tt>S</tt> defined by something like </p>

<pre>&nbsp;&nbsp; SumRandom S = RV1(i1)*r1 - RV2*r2 + RV3(i3)*RV4 + r0;</pre>

<p>has the obvious meaning. <tt>RV1(i1)</tt> means that the sum of <tt>i1</tt> independent
values from <tt>RV1</tt> should be used. Note that <tt>RV1*RV1</tt> means the product of
two independent numbers generated from <tt>RV1</tt>. Remember that <tt>SumRandom</tt> is
slow if the number of terms or copies is large. I support the four arithmetic operators <tt>+</tt>,
<tt>-</tt>, <tt>*</tt> and <tt>/</tt> but cannot calculate the means and variances if you
divide by a random variable. </p>

<p>Use <tt>SumRandom</tt> to quickly set up simple combinations of the existing
generators. But if the combination is going to be used extensively, then it is probably
better to write a new class to do this. </p>

<p>Example: <i>normal</i> with mean = 10, standard deviation = 5: </p>

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

<p>Example: <i>F</i> distribution with <i>m</i> and <i>n</i> degrees of freedom: </p>

<pre>&nbsp;&nbsp; int m, n;
&nbsp;&nbsp; ... put values in m and n
&nbsp;&nbsp; ChiSq Num(m); ChiSq Den(n);
&nbsp;&nbsp; SumRandom F = (double)n/(double)m * Num / Den;
&nbsp;&nbsp; for (int i=0; i&lt;100; i++) cout &lt;&lt; F.Next() &lt;&lt; &quot;\n&quot;;</pre>

<h3><a NAME="mixed"></a><font COLOR="#FF0000">MixedRandom:</font></h3>

<p>This is for mixtures of distributions. Suppose <tt>rv1</tt>, <tt>rv2</tt>, <tt>rv3</tt>
are random number generators and <tt>p1</tt>, <tt>p2</tt>, <tt>p3</tt> are Reals summing
to 1. Then the generator <tt>M</tt> defined by </p>

<pre>&nbsp;&nbsp; MixedRandom M = rv1(p1) + rv2(p2) + rv3(p3);</pre>

<p>produces a random number generator with selects its next random number from <tt>rv1</tt>
with probability <tt>p1</tt>, <tt>rv2</tt> with probability <tt>p2</tt>, <tt>rv3</tt> with
probability <tt>p3</tt>. </p>

<p>Alternatively one can use the constructor </p>

<pre>&nbsp;&nbsp; MixedRandom M(n, prob, rv);</pre>

<p>where <tt>n</tt> is the number of distributions in the mixture, <tt>prob</tt> the Real
array of probabilities, <tt>rv</tt> an array of pointers to random variables. </p>

<p>Normal with outliers: </p>

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

<p>or: </p>

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

<h3><a NAME="permutation"></a><font COLOR="#FF0000">RandomPermutation:</font></h3>

<p>To draw <tt>M</tt> numbers without replacement from <tt>start, start+1, ..., start+N-1</tt>
use </p>

<pre>&nbsp;&nbsp; RandomPermutation RP;
&nbsp;&nbsp; RP.Next(N, M, p, start);</pre>

<p>where <tt>p</tt> is an <tt>int*</tt> pointing to an array of length <tt>M</tt> or
longer. Results are returned to that array. </p>

<pre>&nbsp;&nbsp; RP.Next(N, p, start);</pre>

<p>assumes <tt>M = N</tt>. The parameter, <tt>start<font FACE="Times New Roman">,</font></tt>
has a default value of 0. </p>

<p>The method is rather inefficient if <tt>N</tt> is very large and <tt>M</tt> is much
smaller. </p>

<h3><a NAME="combination"></a><font COLOR="#FF0000">RandomCombination:</font></h3>

<p>To draw <tt>M</tt> numbers without replacement from <tt>start, start+1, ..., start+N-1</tt>
and then sort use </p>

<pre>&nbsp;&nbsp; RandomCombination RC;
&nbsp;&nbsp; RC.Next(N, M, p, start);</pre>

<p>where <tt>p</tt> is an <tt>int*</tt> pointing to an array of length <tt>M</tt> or
longer. Results are returned to that array. </p>

<pre>&nbsp;&nbsp; RC.Next(N, p, start);</pre>

<p>assumes <tt>M = N</tt>. The parameter, <tt>start<font FACE="Times New Roman">,</font></tt>
has a default value of 0. </p>

<p>The method is rather inefficient if <tt>N</tt> is large. A better approach for large <tt>N</tt>
would be to generate the sorted combination directly. This would also provide a better way
of doing permutations with large <tt>N</tt>, small <tt>M</tt>. </p>

<h3><a NAME="extreal"></a><font COLOR="#FF0000">ExtReal</font></h3>

<p>A class consisting of a Real and an enumeration, <tt>EXT_REAL_CODE</tt>, taking the
following values: 

<ul>
  <li>Finite</li>
  <li>PlusInfinity</li>
  <li>MinusInfinity</li>