%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---------------------------------------------------------------------------
% Generating Pseudo-random Numbers ******************************************
% ---------------------------------------------------------------------------
% Submitted:      Dr. Eric Suess,Statistics 6501,Winter 2002
% Submitted Date: January 7st, 2002
% Accepted Date:  January 9th, 2002
% ---------------------------------------------------------------------------
% This is a first draft of the handout "Simulating Discrete Random Variables"
% in LaTeX.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[12]{report}

\usepackage{latexsym}
\usepackage[centertags]{amsmath}
\usepackage{amsfonts}
\usepackage{newlfont}
%\usepackage{graphicx,fullpage,subfigure,pstricks}

\setlength{\oddsidemargin}{-0.4cm}

\setlength{\evensidemargin}{0cm}

\setlength{\topmargin}{-0.75cm}

\setlength{\textheight}{8.9in}

\setlength{\textwidth}{6.5in}

\renewcommand{\baselinestretch}{1}

\setlength{\parskip}{3mm}

\setlength{\parindent}{5mm} % indented paragraphs just don't

%\pagestyle{plain}

\begin{document}
\pagestyle{myheadings} \markright{Statistics 6501, Winter 2002}

\large
\begin{center}
\textbf{Simulating Discrete Random Variables}
\end{center}

\normalsize
\noindent
\textbf{Binomial} \\

One way to simulate a binomial random variable is to
simulate the events of which it is composed. If this binomial
random variable has parameters $n$ and $p$, then we can
independently simulate $n$ events, each with probability $p$ of
being a $success$. The number of $successes$, then, would be the
outcome of the binomial random variable.  The code listed below
uses this method to create an integer valued function that returns
the outcome of the random variable: (Note that in this function
and in all other code in this discussion the type
\texttt{extended} is used.  This is a floating point variable type
that is defined in \textrm{Turbo Pascal}. It has about 20
significant digits and can range from $3.4 \times 10^{-4932}$ to
$1.1 \times 10^{4932}$. This provides for greater accuracy in some
of these applications; this is particularly necessary in the
Poisson function on the next page. I recommend that you take
advantage of any high precision real types that your compiler
supports.  In \textrm{Fortran} and \textrm{C} the type \texttt{double}
carries about 16 significant digits.)

\begin{verbatim}
function binomial(n : integer; p : double): integer;

var i : integer;
    successes : integer;

begin
successes := 0;
for i := 1 to n do
   if random < p then inc(successes);
binomial := successes;
end;
\end{verbatim}

\noindent
\textbf{Geometric} \\

Here again we will use the method of simulating the
events that make up the random variable.  For a geometric random
variable with parameter $p$ we simulate independent events, each
with probability $p$ of a $success$ occurring, until we observe a
$success$. The number of tries it takes to get the first $success$
is the outcome of the random variable.

\begin{verbatim}
function geometric(p : double): integer;

var i : integer;
    success : boolean;  (* becomes true if a success occurs *)

begin
success := false;
i := 0;
while not success do
  begin
  success := random < p;
  i := i + 1;
  end;
geometric := i;
end;
\end{verbatim}

\noindent
\textbf{Poisson} \\

The listing below is a function that will return the
outcome of a poisson random variable with parameter $\lambda$.  We
will not be able to discuss the method used here until later in the
course.

\begin{verbatim}
function poisson(lambda : double): longint;

var i : longint;
    product : double;
    compare : double;

begin
compare := exp(-lambda);
product := random;
i := 0;
while product > compare do
  begin
  product := product * random;
  i := i + 1;
  end;
poisson := i;
end;
\end{verbatim}

\newpage        % starts new page

\noindent
\textbf{Application} \\

The code listed below generates 1000 observations from a
binomial random variable with $n = 20$ and $p = 0.10$.  It keeps
track of the min and max of all the outcomes.  Also it keeps track
of $\sum_{i=1}^{1000}x_i$ and $\sum_{i=1}^{1000}{x_i^2}$ and
uses these to calculate the sample mean (usually denoted by
$\bar{x}$) and the sample variance (usually denoted by $s^2$).
These have the following computational formulas:

\large
\begin{center}
\begin{tabular}{cc}
$\bar{x} = \frac{\sum_{i=1}^{m} x_i}{m}$ & \hspace{1.5in}
$s^2 = \frac{\sum_{i=1}^{m}{x_i^2} - \frac{{(\sum_{i=1}^{m}{x_i})}^2}{m}}{m-1}$
\end{tabular}
\end{center}
\normalsize

where here $m = 1000$.  If $m$ is large (1000 is fairly large)
then these shouldn't be too far from the true mean and true
variance ($\mu = np$ and $\sigma^2 = np(1-p)$).

\begin{verbatim}
program GenerateDiscreteDistributions;

var x : double;
    Min, Max : integer;
    SampleMean, SampleVariance : double;
    Sum, SumSq : double;
    i : integer;

(* insert random function here *)

(* insert binomial function here *)

begin
Randomize; (* Seeds Turbo Pascal's random number generator *)
Min := 50;
Max := 0;
Sum := 0;
SumSq := 0;
For i := 1 to 1000 do
begin
  x := Binomial(20,0.10);
  Sum := Sum + x;
  SumSq := SumSq + x * x;
  if x < Min then Min := x;
  if x > Max then Max := x;
end;
SampleMean := Sum/1000;
SampleVariance := (SumSq - (Sum * Sum)/1000)/1000;
WriteIn('Sample mean = ', SampleMean:1:6);
WriteIn('Sample variance = ', SampleVariance:1:6);
WriteIn('min = ', Min:1:0);
WriteIn('max = ', Max:1:0);
\end{verbatim}

\end{document}
% ------------------------------------------------------------------------
