Applied
Probability and Stochastic Processes
Lecture Notes for 577/578 Class

Wlodzimierz Bryc
Department of Mathematics
University of Cincinnati
P. O. Box 210025
Cincinnati, OH 45221-0025
E-mail: brycw@math.uc.edu

Created: October 22, 1995
Modified: 1996
Printed: Mar 13, 1998
© 1995 Wodzimierz Bryc
WWW version

Course description

This course is aimed at students in applied fields and assumes a prerequisite of calculus. The goal is a working knowledge of the concepts and uses of modern probability theory. A significant part of such a ``working knowledge" in modern applications of mathematics is computer-dependent.

The course contains mathematical problems and computer exercises. Students will be expected to write and execute programs pertinent to the material of the course. No programming experience is necessary or assumed. But the willingness to accept a computer as a tool is a requirement.

For novices, BASIC programming language, (QBASIC in DOS, Visual Basic in Windows, or BASIC on Macintosh) is recommended. BASIC is perhaps the easiest programming language to learn, and the first programming language is always the hardest to pick.

Programs in QBASIC 4.5 on IBM-compatible PC and, to a lesser extend, programs on Texas Instrument Programmable calculator TI-85, and Windows_TM programs in Microsoft Visual Basic 3.0 are supported. This means that I will attempt to answer technical questions and provide examples. Other programming languages (SAS, C, C++, Fortran, Cobol, Assembler, Mathematica, LISP, T_EX (!), Excel, etc.) can be used, but I will not be able to help with the technical issues.

Contents of the course (subject to change)

• [577] Basic elements of probability. Poisson, geometric, binomial, normal, exponential distributions. Simulations. Conditioning, characterizations.

  Moment generating functions, limit theorems, characteristic functions. Stochastic processes: random walks, Markov sequences, the Poisson process.

• [578] Time dependent and stochastic processes: Markov processes, branching processes. Modeling

  Multivariate normal distribution. Gaussian processes, white noise. Conditional expectations. Fourier expansions, time series.

Supporting materials

This text is available through the Internet in PDF, PostScript, or DVI formats. A number of other math related resources can be found on WWW. Also available are supporting BASIC program files.

Auxiliary textbooks:

• W. Feller, An Introduction to Probability Theory, Vol. I, Wiley 1967. Vol II, Wiley, New York 1966
  Volume I is an excellent introduction to elementary and not-that-elementary probability. Volume II is advanced.

• W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical recipes in C, Cambridge University Press, New York
  A reference for numerical methods: C-language version.
• J. C. Sprott Numerical recipes Routines and Examples in BASIC, Cambridge University Press, New York 1992
  A reference for numerical methods: routines in QuickBasic 4.5 version.

• H. M. Taylor & S. Karlin, An introduction to stochastic modeling, Acad. Press, Boston 1984
  Markov chains with many examples/models, Branching processes, Queueing systems.

• L. Breiman, Probability and Stochastic Processes: with a view towards applications, Houghton Mifflin, Boston 1969.
  includes Markov chains and spectral theory of stationary processes.

• R. E. Barlow & F. Proschan, Mathematical Theory of Reliability, SIAM series in applied math, Wiley, New York 1965.
  Advanced compendium of reliability theory methods (mid-sixties).

• S. Biswas, Applied Stochastic Processes, Wiley, New York 1995
  Techniques of interest in population dynamics, epidemiology, clinical drug trials, fertility and mortality analysis.

• J. Higgins & S. Keller-McNulty, Concepts in Probability and Stochastic Modeling Duxbury Press 1995
  Covers traditional material; computer simulations complement theory.

• T. Harris, The Theory of Branching Processes Reprinted: Dover, 1989.
  A classic on Branching processes.
• H. C. Tjims, Stochastic models. An algorithmic approach, Wiley, Chichester, 1994.
  Renewal processes, reliability, inventory models, queuieing models.

Conventions

Exercises

Programming

Copyright and License for 1996 version

This textbook is copyrighted © W. Bryc 1996. All rights reserved. GTDR.

Academic users are granted the right to make copies for their personal use, and for distribution to their students in academic year 1995/1996.

Reproduction in any form by commercial publishers, including textbook publishers, requires written permission from the copyright holder.

WWW version© 1996 W. Bryc

Part 1
Probability

Chapter 1
Random phenomena

Definition of a Tree: A tree is a woody plant with erect perennial trunk of at least 3.5 inches (7.5 centimeters) in diameter at breast height (4¹/₂ feet or 1.3 meters), a definitely formed crown of foliage, and a height of at least 13 feet (4 meters).
The Auborn Society Field Guide to North American Trees.

This chapter introduces fundamental concepts of probability theory; events, and their chances. For the readers who are familiar with elementary probability it may be refreshing to see the computer used for counting elementary events, and randomization used to solve a deterministic optimization problem.

The questions are

• What is ``probability"?
• How do we evaluate probabilities in real-life situations?
• What is the computer good for?

1.1  Mathematical models, and stochastic models

To begin with, we start with a truism. Real world is complicated, often to a larger degree than scientists will readily admit. Most real phenomena have multi-aspect form, and can be approached in multiple ways. Various questions can be asked. Even seemingly the same question can be answered on many incompatible levels. For instance, the generic question about dinosaur extinction has the following variants.

• Why did Dino the dinosaur die? Was she bitten by a poisonous rat-like creature? Hit by a meteorite? Froze to death?
• Why did the crocodiles survive to our times, and tyranosaurus rex didn't?
• What was the cause of the dinosaur extinction?
• Was the dinosaur extinction an accident, or did it have to happen? (This way, or the other).
• Do all species die out?

Theses questions are ordered from the most individual level to the most abstract. The reader should be aware that probability theory, and stochastic modelling deal only with the most abstract levels of the question. Thus, a stochastic model may perhaps shed some light whether dinosaurs had to go extinct, or whether mammals will go extinct, but it wouldn't go into details of which comet had to be responsible for dinosaurs, or which is the one that will be responsible for the extinction of mammals.

It isn't our contention that individual questions have no merit. They do, and perhaps they are as important as the general theories. For example, a detective investigating the cause of a mysterious death of a young woman, will have little interest in the ``abstract statistical fact" that all humans eventually die anyhow. But individual questions are as many as the trees in the forest, and we don't want to overlook the forest, either.

Probabilistic models deal with general laws, not individual histories. Their predictions are on the same level, too. To come back to our motivating example, even if a stochastic model did predict the extinction of dinosaurs (eventually), it would not say that it had to happen at the time when it actually happened, some 60 million years ago. And the more concrete a question is posed, say if we want to know when Dino the dinosaur died, the less can be extracted from the stochastic model.

On the other hand, many concepts of modern science are define in statistical, or probabilistic sense.(If you think this isn't true, ask yourself how many trees do make a forest.) Such concepts are best studied from the probabilistic perspective. An extreme view is to consider everything random, deterministic models being just approximations that work for small levels of noise.

1.2  Events and their chances

Let \cal M be a s-field of its subsets, called the events. Events A,B Î \cal M model sentences about the outcomes of the experiment to which we want to assign probabilities. Under this interpretation, the union AÈB of events corresponds to the alternative, the intersection AÇB corresponds to the conjunction of sentences, and the complement A¢ corresponds to the negation of a sentence. For A,B Î \cal M, A\B: = AÇB¢ denotes the set-theoretical difference.

For an event A Î \cal M the probability Pr(A) is a number interpreted as the degree of certainty (in unique experiments), or asymptotic frequency of A (in repeated experiments). Probability Pr(A) is assigned to all events A Î \cal M, but it must satisfy certain requirements (axioms). A set function Pr(·) is a probability measure on (W,\cal M), if it fulfills the following conditions:

1. 0 £ Pr(A) £ 1
2. Pr(W) = 1
3. For disjoint¹ Pr(AÈB) = Pr(A)+Pr(B).
4. If A_n are such that Ç_{n ³ 1} A_n = Æ and A₁ É A₂ É ... É A_n É A_n+1 É ... are decreasing events, then

   Pr (An)® 0.

   (1)

Probability axioms do not determine the probabilities in a unique way. The axioms provide only minimal consistency requirements, which are satisfied by many different models.

1.2.1  Uniform Probability

The uniform assignment of probability involves counting. For small sample spaces this can be accomplished by examining all cases. Moderate size sample spaces can be inspected by a computer program. Counting arbitrary large spaces is the domain of combinatorics. It involves combinations, permutations, generating functions, combinatorial identities, etc. Short recalls the most elementary counting techniques. Problem 1 Three identical dice are tossed. What is the probability of two of a kind?

The following BASIC program inspects all outcomes when five dice are rolled, and counts how many are ``four of a kind".

PROGRAM yahtzee.bas
'

'declare function and  variables
DECLARE FUNCTION CountEq! (a!, b!, c!, d!, e!)

'prepare screen
CLS
PRINT "Listing four of a kind outcomes in Yahtzee..."

'*** go through all cases
FOR a = 1 TO 6: FOR b = 1 TO 6: FOR c = 1 TO 6: FOR d = 1 TO 6: FOR e = 1 TO 6
   IF CountEq(a + 0, b + 0, c + 0, d + 0, e + 0) = 4 THEN
       PRINT a; b; c; d; e; : ct = ct + 1
       IF ct MOD 5 = 0 THEN PRINT :  ELSE PRINT "|";
   END IF
NEXT: NEXT: NEXT: NEXT: NEXT

'print result
PRINT
PRINT "Total of "; ct; " four of a kind."

FUNCTION CountEq (a, b, c, d, e)
'*** count how many of five numbers are the same
DIM x(5)
x(1) = a
x(2) = b
x(3) = c
x(4) = d
x(5) = e
max = 0
FOR j = 1 TO 5
    ck = 0
    FOR k = 1 TO 5
        IF x(j) = x(k) THEN ck = ck + 1
    NEXT k
    ck = -ck
    IF ck > max THEN max = ck
NEXT j
'assign value to function
CountEq = max
'

Here is a portion of its output:

The program is written explicitly for tossing five dice, and you may want to modify it to answer similar question for any number of dice, and an arbitrary k-of-a-kind question.

[ 1 Run and time YAHTZEE.BAS. Then estimate how long a similar problem would run if the question involved tossing 15 fair dice. The answer depends on your computer, and the software. Both Pascal and C-programs seem to run on my computer about 15 times faster than the (compiled) QuickBasic. ANS: Running time would take years!

Exercise 1.2.1 shows the power of old-fashioned pencil-and-paper calculation. Problem 2 Continuing Problem 1.2.1, suppose now n identical dice are tossed. What is the probability of n-1 of a kind? ANS: [5 n/( 6^n-1)].

1.2.2  Geometric Probability

For bounded subsets W Ì IR^d, put

Geometric probability usually involves multivariate integrals. Example 1 A point is selected from the 32 cm wide circular dartboard. The probability that it lies within the 8cm wide bullseye is [1/ 16].

Example 2 A needle of length 2l < 2 is thrown onto a paper ruled every 2 inches. The probability that the needle crosses a line is 2l/p.

Analysis of Example 1.2.2 is available on WWW.

[ 2 Test by experiment (or simulation on the computer) if the above two examples give correct answers. (Before writing a program, you may want to read Section first.)

Example 3 Two drivers arrive at an intersection between 8:00 and 8:01 every day. If they arrive within 15 seconds of each other, both cars have to stop at the stop-sign. How often do the drivers pass through the intersection without stopping?

[ 1 Continuing Example 1.2.2: What if there are three cars in this neighborhood? Four? How does the probability change with the number of cars? At what number of users a stop light should be installed?

1.3  Elementary probability models

1.3.1  Consequences of axioms

1.3.2  General discrete sample space

For a finite or countable set W Ì IN and a given summable sequence of non-negative numbers a_n put

p_k = [(a_k)/( å_{n Î W} a_n)] and rewrite (8) as Pr(A) = å_{n Î A} p_n.

Formula (8) generalizes the uniform assignment (2), which corresponds to the choice of equal weights a_k = 1. At the same time it is more flexible² and applies also to infinite sample spaces. Example 4 Let W = IN and put p_k = [1/( 2^k)]. Then the probability that an odd integer is selected is Pr(Odd) = å_{j = 0}^¥ 2^-2j-1 = ²/₃. ()

Table list the most frequently encountered discrete probability assignments.

Name	W	Probabilities pk
Binomial	{0,...,n}	pk = (nk)pk(1-p)n-k
Poisson	ZZ\scriptscriptstyle +	pk = e-l[(lk)/ k!]
Geometric	IN	pk = p(1-p)k-1
Equally likely outcomes	{x1,...,xk}	pk = 1/k

Problem 3 For each of the choices of numbers p_k in Table 1.1 verify that (8) indeed defines a probability measure.

The reasons behind the particular choices of the expressions for p_k in Table 1.1 involve modeling.

1.3.3  General continuous sample space

1.4  Simulating discrete events

The first step in a simulation is to decide how to generate ``randomness" within a deterministic algorithm in a computer program. Programming languages, and even calculators, provide a method for generating consecutive numbers from the interval (0,1) ``at random". For instance, BASIC instruction³ PRINT RND(1) returns different number at every use. We shall assume that the reader has access to a of generating uniform numbers from the unit interval (0,1). These are called pseudo-random numbers, since the program usually begin the same ``random" sequence at every run, unless special precautions⁴ are taken.

Once a pseudo-random number from the interval (0,1) is selected, an event that occurs with some known probability p can be simulated by verifying if { RAND(1) < p} occurs in the program. For instance, the number of heads in a toss of a 1,000 fair coins is simulated by the following BASIC program.

PROGRAM tosscoin.bas
'

'Simulating 1000 tosses of a fair coins

H = 0
FOR n = 1 TO 1000  'main loop
        IF RND(1) < .5 THEN H = H + 1
NEXT n
'print final message
PRINT "Got "; H; " heads this time"
END
'

Here is its output:Got 520 heads this time

A simple method for simulating an outcome on a six-face die is to take the integer part of a re-scaled uniformly selected number INT(1+6*RND(1)). Can you use this to write a simulation of a roll of 5 dice? Such simulations are often used to evaluate probabilities empirically as a substitute for a real empirical study.

[ 3 Write the simulation (as opposed to deterministic inspection of all sample points on page pageref) to estimate how often the event ``four of a kind" occurs in a roll of five dice.

[ 2 Run the (modified) coin tossing program in a loop and to answer the following questions.

• In a 100 tosses of a coin, how often does less than 55 heads occur?
• In a 100 tosses of a coin, how often does less than 80 heads occur?
• Can you sketch the curve⁵ that represents the probability of less than x coins in n = 100 tosses of a fair coin?

More complicated objects are often of interest in simulations. For instance we may want to draw two cards from the deck of 52. One possible way to do it is to number the cards, and select two numbers a,b.

1. Select the first card a=INT(RND(1)*52+1) as a random integer between 1 and 52.
2. Select the second card b=INT(RND(1)*52+1) in the same way.
3. Compare a,b

   1. If a = b then repeat step 2
   2. Otherwise, got two different cards a ¹ b at random

1.4.1  Generic simulation template

In more advanced exercises you may want to estimate probabilities that aren't known. In such cases it is always a good idea to run simulations of various lengths and compare the results. In this section we briefly discuss how such a simulation can be organized in a way that promotes multiple uses of the same program.

The key is organizing the programs carefully into manageable blocks of small size. Modern BASIc is a structural programming language. The generic program to study the effects of the length of simulation on its output can be written as follows

' PROGRAM Generic.bas
'Generic Simulator
'Size is simulation size varied from Min=100 to Max=10000
For Size 100 to 10000 Step 100
Simulate(Size, Result)
Print "Simulation size="; Size ; Öutput="; Result
Next Size
End

The actual simulation is performed by

SUB Simulate (SizeRequested, Result)
'Runs requested number of simulations and returns average
'Trial numbers consecutive simulations from 1 to SizeRequested
For Trial=1 to SizeRequested
SimulateOne(Score)
Result=Result+Score
Next Size
'Most simulations return averages of single trials
Result=Score/Size
End SUB

The actual modeling is performed in another SUB, which in the generic program we named SimulateOne. This SUB may be as simple as simulating a toss of a single coin

SUB SimulateONe(Outcome)
'simulate One occurrence, return numerical outcome
OutCome=0
if RND(1)<1/2 THEN Outcome=1
END SUB

Or it can be as complicated as we wish. The example below simulates a toss of five dice, and uses previously introduced function CountEq(a,b,c,d,e). four-of-a-kind.

SUB SimulateONe(Outcome)
'simulate One occurrence, return numerical outcome
d1=int(RND(1)*6+1)
d2=int(RND(1)*6+1)
d3=int(RND(1)*6+1)
d4=int(RND(1)*6+1)
d5=int(RND(1)*6+1)
IF CountEq(d1,d2,d3,d4,d5)=4 THEN Outcome=1
END SUB

1.4.2  Blind search

[ 3 Write a blind-search program to find the maximum of a function.

• Organize your program so that the function can easily be changed - but for now use the one you are quite familiar with, like 100-(x-300)², or 300e^-(x-30)2sin(200x).

• If you are looking for more challenge, do the same for three variables. Write a blind-search program to find a maximum of a function like 300e^-(x-30)2sin(200x+400y-z)+400e^-(y-70)2cos(400x+200y+z) over the ball x²+y²+z² £ 1000.
• As a further complication, try to find a maximum of a function that has two local maxima, and the region isn't convex. (This is an almost hopeless task for gradient methods!)

1.4.3  Application: Traveling salesman problem

Planning such a tour is easy by hand for 3-4 cities. For longer tours some help is needed. To check the performance of the blind search, you may want to know what the usual algorithms involve. A greedy method starts with the shortest distance, and then keeps going to the closest city not yet visited. Another heuristic method is to to select a pair of closest cities and accept the best (shortest) connection among those that do not complete a shorter loop, and do not introduce a spurious third connection to a city. Eventually, the disjoint pieces will form a path that often can be further improved upon inspection.

The program is longer but not at all sophisticated - it just selects paths at random. Notice that a solution to Exercise 1.4 - how to generate a random permutation - is given in one of the subprograms (which one?). The method for the latter is simple-minded - the algorithm attempts to place consecutive numbers at random spots until an empty spot is selected.

The following is the main part of the program. You can use it as a template in designing your own version of Blind search programs. The full code with SUBs is online in RANDTOUR.BAS.

'****
CLS
'**** get number of cities (no choice which) from user
LOCATE 2, 1
INPUT ; "Shortest distance between how many cities?", n

'*** initialize program
CLS
nMax = 19 ' current data size. Make sure not exceeded!
IF n > nMax THEN n = nMax

'declare arrays
DIM SHARED dist(nMax, nMax) ' matrix of distances
DIM SHARED city(nMax) AS STRING
DIM P(n), BestP(n)

'read distances
CALL AssignDistances(nMax)

'initial permutation
FOR j = 1 TO n
        P(j) = j
        BestP(j) = j
NEXT j
'initial length of trip
MinLen = PathLen(P())

'*** main loop
DO  'infinite loop till user stops
'count trials
no = no + 1
        '*** interacting with user
        'check if user pressed key to stop
        k$ = INKEY$
        IF k$ > "" THEN EXIT DO  'exit infinite loop
        'display currrent progress
        display (no)
        '*** get any path
        CALL GetPermutation(P())
        x = PathLen(P())
        IF x < MinLen THEN
                'Better path found, so memorize and display
                Dlen = MinLen - x
                MinLen = x
                '*Memorize best order and print 
                FOR j = 1 TO n
                        BestP(j) = P(j)
                        PRINT city(P(j)); "->";
                NEXT j
                'Finish printing
                PRINT city(P(1))
                PRINT "Best so far: "; MinLen
                PRINT "Progress rate "; Dlen / (no - Slen); " miles per trial"
                Slen = no
        END IF
LOOP
'Print final message
CLS
PRINT "ALPHABETIC TOUR OF FIRST "; n; " CITIES in the USA"
PRINT "Blind Search Recommended Route found in "; Slen; "-th search"
FOR j = 1 TO n - 1
        PRINT city(BestP(j)); "-->";
NEXT j
PRINT city(BestP(n)); "-->"; city(BestP(1))
PRINT "Total distance: "; MinLen
LOCATE 22, 40
PRINT "(Distances subject to change)"
END

When you run this program on larger sets of cities, you will notice that the program is not fast. One possible improvement in the design of this program is to modify the randomization to be less likely to pick long paths. For instance, you can attempt to modify paths that are known to be short, or weight the modifications by lengths of resulting paths. Such methods are actually in use in image restoration problems (simulated annealing), see page .

1.4.4  Improving blind search

• Blind search programs are easy to write, if you know how to code the main function to randomize.
• Blind search always gives ``answers"
• It is difficult to judge how good an answer is.
• In situations that we do know the answer, blind search takes long time to reach it.

It is possible to improve on the last aspect without complicating the program much. The idea is to make random modifications of the currently best found value. For example, in a one-dimensional maximization of a function f(x), we would do the following steps

1. Pick an initial ``best-so-far" point x₀ and compute initial value y₀ = f(x₀)
2. Select at random x₁ in the ``neighborhood" of x₀ and compute y₁ = f(x₁)
3. Compare y₀, y₁.

   1. If y₁ < y₂ then repeat Step 2.
   2. If y₁ ³ y₀ then make x₁ the new ``best-so-far" y₀: = y₁,x₀: = x₁. Then repeat Step 2.

4. Stop the program at user request, or when no changes to y₀ occur for prolonged number of attempts to improve it.

We want to allow for the chance of checking points far away from the ``best-so-far" answer. But we don't want this to happen too often, because x<sub>0</sub> might be rather close to the optimum. The tradeoffs are that the program will tend to follow ``direct path" to the maximum, but the danger is that it will get stuck longer in local maxima.

Improved blind search with time/state dependent randomization is actually implemented within RANDTOUR. It is commented out in the version on the disk, so that it isn't operational. To make it active, uncomment the call to ImproveBest as a replacement for GetPermutation.

1.4.5  Random permutations

Program RANDTOUR.BAS selects permutations at random only in its ``simplest" variant. Here are a few examples of problems that require selecting random permutations.

• Card games:
  • Poker hand: Select 5 cards at random from a deck of cards.
  • Poker (2 players): Select 10 cards at random from a deck of cards.
  • Bridge: Split 52 cards into four groups

• Analyzing statistical experiments:

  Suppose there are 7 items hidden under 12 cups, and a person is allowed to try to find them. How often all seven will be recovered by pure lack (as opposed to, say, parapsychic abilities?

SUB GetPermutationFast(P())
'Put random integers into P()
n=UBOUND(P) 'how many entries
'this loop can be omitted is we are sure that P(j) list all the numbers we want
FOR j=1 TO n
P(j)=j
next j
for j=1 to UBOUND(P) 'not n as n changes in the loop
k=INT(RND(1)*n+1)

SWAP P(n), P(k)
n=n-1
next j

[ 5 There are many incompatible measures of ``quality" of an algorithm. One of the ``objective" criteria is the number of ``If" verifications. Another ``objective criterion" might be the number of calls to a function. A ``subjective" criterion, which depends on the hardware and circumstances, is timing.

Does SUB GetPermutationFast deserve adverb Fast in its name?

[ 4 The Subset-Sum Problem is stated as follows:

Let S be a set of positive integers and let t be a positive integer. Decide of there is a subset S¢ Ì S such that the sum of integers in S¢ is exactly t.

The associated optimization problem is to find a subset S¢ Ì S whose elements' sum is the largest but doesn't exceed t. This optimization problem is NP-complete, ie it isn't known if there is a polynomial time algorithm (polynomial in the size of S) to find S¢.

Investigate how the blind search will do on sets S selected at random, and on sets S constructed in more regular fashion like arithmetic progression S = {a, 2a, 3a, ...}, geometric progression S = a,a², a³,....

1.5  Conditional probability

To formalize this idea, suppose B is an event such that Pr(B) ¹ 0. The last condition merely says that B is an event that does have some chance of occurring. Conditional probability of event A given event B is denoted by Pr(A|B). It is defined as

Conditional probability satisfies the axioms of probability, and Pr(A|B) = 0 if A,B are disjoint. In particular, Pr(A¢|B) = 1-Pr(A|B), Pr(B|B) = 1.

The easiest way to find Pr(A|B) by simulations is to repeatedly simulate the complete experiment, discarding all the outcomes except the ones resulting in B.

[ 6 Suppose we toss repeatedly a fair coin, and the ``success" is to get heads.

Use computer simulations to find the conditional probability that the very first trial was successful, if 10 consecutive (and independent) trials resulted in 8 successes.

Try to answer the same question under the condition that 50 independent trials resulted in 40 successes.

You should notice that it takes forever to simulate events that happen rarely. Section indicates one possible way out of this difficulty.

1.5.1  Properties of conditional probability

Example 5 Suppose we have a deck of 52 cards numbered 1 through 52. Since it isn't obvious how to simulate selecting 5 cards without replacement, we may want to select them with replacement instead. Let A denote the event that all five cards are different. What is the probability of A?

We may perform the experiment sequentially, drawing one card at a time. Let A_k denote the event that the k consecutive draws resulted in different cards. Then A = A₅ Ì A₄ Ì ... Ì A₁. Moreover, Pr(A₁) = 1.

Clearly, Pr(A₂|A₁) = [51/ 52], so Pr(A₂) = Pr(A₂ÇA₁) = Pr(A₂|A₁)Pr(A₁) = [51/ 52]. Similarly, Pr(A₃) = Pr(A₃ÇA₂) = Pr(A₃|A₂)Pr(A₂) = [50/ 52][51/ 52]. Continuing this we get Pr(A₅) = [51 50 49 48/( 52⁴)] » 0.82.

The following identities are also of interest.

1. Path Probability: Pr(Ç_{k = 1}ⁿ A_k) = Õ_{k = 1}ⁿPr(A_k|Ç_{j = 1}^k-1 A_j)
2. Bayes theorem: Pr(A|B) = [(Pr(B|A)Pr(A))/Pr(B)]
3. Total probability formula: If {B_n} are pairwise disjoint and exhaustive, ie. Pr(A_iÇB_j) = Æ for i ¹ j and ÈB_n = W, then

   Pr (A) = å n  P(A|Bn) Pr (Bn).

   (10)

Example 6 A lake has 200 white fish and a 100 black fish, and a nearby pond contains 20 black fish and 10 white ones. No other fish live there.

A fish is selected at random from the lake and moved to the pond. Then a fish is selected from the pond and moved back to the lake. What is the probability that all fish in the pond are black?

1.5.2  Sequential experiments

Denoting by \cal P the generic path A₁ÇA₂Ç...ÇA_n, and by \cal P(k) = A_k we have the following path integral formula for the probability of an event \cal F specifying the outcome of the complete experiment, and consisting of many paths \cal P.

Example 7 Suppose that the double transfer operation from Example 1.5.1 was repeated twice. That is, a random selection was done four times. What is the probability that all fish in the pond are black?

[ 8 Check by simulation how the proportion of black fish in the pond changes when the random transfers from Example 1.5.1 are performed repeatedly for a long time.

1.6  Independent events

Two events A,B are independent, if the conditional probability is the same as unconditional, Pr(A|B) = Pr(A). This is stated in multiplicative form which exhibits symmetry and includes trivial events⁸ Definition 1 Events A, B are independent if Pr(AÇB) = Pr(A)Pr(B).

Independence captures the intuition of non-interaction, and lack of information. In modeling it is often assumed rather than verified. For instance, we shall assume that the events generated by consecutive outputs of the random generator are independent. We also assume that tosses of a coin (fair, or not!) are independent.

Beginners sometimes confuse disjoint versus independent events. Exclusive (ie. disjoint) events capture the intuition of non-compatible outcomes. Not compatible outcomes cannot happen at the same time. This is not the same as independent outcomes. If A, B are disjoint and you know that A occurred, then you do know a lot about B. Namely you know that B cannot occur. Thus there is an interaction between A and B. Knowing whether A occurred influences chances of B, which is not possible under independence.

Independence (or, more properly, mutual stochastic independence) of families of events is defined by requesting a much larger number of multiplicative conditions. The reason behind is Theorem , which provides a very convenient tool. Definition 2 Events A₁, A₂, ..., A_n are independent, if Pr(Ç_{j Î J} A_j) = Õ_{j Î J} Pr(A_j) for all finite subsets J Ì IN.

Example 8 A coin is tossed repeatedly. Find the probability that heads appears for the first time on the fourth toss.

Problem 4 SUB GetPermutation from the program RANDTOUR.BAS selects numbers between 1 and n at random until it finds a number not yet on the list. Then it ads the number to its list, and repeats the process.

1. What is the probability that the second number added to the list required more than k attempts?
2. What is the probability that the last number added to the list required more than k attempts?

Another important concept is the conditional independence. For example, many events in the past and in the future are dependent. But in many mathematical models, past and future are independent conditionally on the present situation. In such a model future depends on past only through present events!

Definition 3 Let C be a non-trivial event. Events A,B are C-conditionally independent if Pr(AÇB|C) = Pr(A|C)Pr(B|C).

1.6.1  Random variables

Random variables are introduced for convenient description of experiments with numerical outcomes. (The other option is to select W Ì IR, or W Ì IR^d.) If we want to run computer simulations, we need to represent even non-numerical experiments (like tossing coins) in numerical terms anyhow. Thus the language of random variables becomes the natural extension of elementary probability theory, expressing many of the same concepts in a little different language.

A random variable is the numerical quantity assigned to every outcome of the experiment. In mathematical terms, random variable is a function X:W® IR with the property that sets {w Î W: X(w) £ a} are events in \cal M for all a Î IR. Recall that the last conditions means that we may talk about probabilities of events {w Î W: X(w) £ a}.

Probabilities for a one-dimensional r. v. are determined by the cumulative distribution function

Cumulative distribution function can be used to express probabilities of intervals Pr(a < X £ b) = F(b)-F(a). Since probability is continuous, (1) we can also compute Pr(X = a) = lim_{b® a⁺}Pr(a < X £ b) = F(a⁺)-F(a). The right hand side limit F(a⁺) exists, as F is a non-decreasing function.

Example 9 Suppose F(x) = (1-e^-x)Ú0. Then Pr(|X-2| < 1) = e^-1-e^-3.

In probability theory we are concerned with probabilities. Random variables that have the same probabilities are therefore considered equivalent. We write X @ Y to denote the equality of distributions, ie. Pr(X Î U) = Pr(Y Î U) for all Borel sets U Ì IR (say, all intervals U).

Vector valued r. v. are measurable the functions W® IR^d. In the vector case we also refer to X = (X₁,... X_d) as the d-variate, or multivariate, random variable.

We will use the ordinary notation for sums and inequalities between random variables. There is however a word of caution. In probability theory, equalities and inequalities between random variables are interpreted almost surely. For instance X £ Y+1 means Pr(X £ Y+1) = 1; the latter is a shortcut that we use for the expression Pr({w Î W: X(w) £ Y(w)+1}) = 1. Problem 5 Show that F(x) = Pr(X £ x) is right continuous: lim_x®a⁺F(x) = F(a).

1.6.2  Binomial trials

A binomial experiment, called also binomial trials, consists of the sequence of simpler identical experiments that have two possible outcomes each. The independent events S_j represent successes in consecutive experiments. We assume that we have an infinite sequence of events S₁,S₂,... S_k,... that are independent and have the same probability p = Pr(S_j). We denote by F_j = S_j¢ the failure in the j-th experiment, and put q = 1-p.

Two important random variables are associated with the binomial experiment are the number X of successes in n trials, and the number T of trials until first success. Example 10 The probability that number X of successes in n trials is k is Pr(X = k) = (_kⁿ)p^kq^n-k. (Here k = 0,..., n.)

Example 11 The probability of more than n attempts needed for the first success is Pr(T > n) = qⁿ. The probability that first success occurs at the n-th trial is Pr(T = n) = pq^n-1 (geometric).

Example 12 Geometric distribution has lack of memory property: Pr(T > n+k|T > n) = Pr(T > k).

Random variables are often described solely in terms of cumulative distribution function F(x), or formulas for Pr(X = x) without reference to the underlying probability space W. For instance, the number of minutes T that we spend waiting for a bird to come to the bird feeder at the back of my house is random, and I believe Pr(T = n) = pq^n-1 because Pr(T > n+k|T > n) = Pr(T > k).

It is intuitively obvious that on average we get np successes in n trials. It is perhaps less obvious¹⁰ that on average we need 1/p trials to get the first success.

[ 9 Write a simulation program to verify the claims about the averages for several values of p.

Example 13 The probability that in 2n tosses of a fair coin, half are heads is [(2n)!/( 4ⁿ(n!)²)] » [1/( Ö{ pn})]® 0 as n®¥. The latter isn't easy to prove, but the computer printout is quite convincing, see Table (note that [1/( Ö{p})] » 0.5642).

2n	Pr(X = n)	Frequency in 1000 trials	ÖnPr(X = n)
100	0.07959	0.08200	0.56278
300	0.04603	0.06100	0.56372
500	0.03566	0.03700	0.56391
700	0.03015	0.02200	0.56399

1.7  Further notes about simulations

By now you should have written some simple simulation programs, and printed out the results. It is perhaps a good moment to pause and consider what are the aspects of simulations that we are interested in.

In general, we would like to get answers to questions that we don't know how to answer in any other way. But before we do that, we should develop some intuition on the cases that can check the answers. Therefore we begin with simulation of probabilities or averages that are known.

• [A] Simulation of probabilities/ averages that are known should address the following questions.

  1. How close the simulation answers are to the theoretical answers? Print them side-by-side.
  2. How large the simulation should be? Is it worth to change simulation size from 1,000 to 10,000 trials? In order to answer this question, your simulation has to provide ``relative" rather than absolute answers. (Answers of the form ``got 32 heads" are meaningless as they depend on simulation size!)
  3. How do the answers change as we change the parameters? If you did a simulation of the fair coin, you could change the probability p from the usual value ¹/₂.

• [B] The next natural step is to extend models that we know how to handle both theoretically and by simulations to cover aspects that aren't easily accessible by theory. The sample questions involve

  1. How would the answers change, if we allow perhaps more realistic assumptions in the model? As an example, suppose that we would like to model the birthday problem with people born non-uniformly throughout the year. Which way would you expect the answer to change?
  2. What are typical errors of a simulation of size n? How can we estimate the accuracy of the answer without having the exact answer to compare it to? Chapter gives theoretical basis for such estimates.

1.8  Questions

Problem 7 [Exercise] A die is thrown until an ace turns up.

Assuming the ace doesn't turn up on the first throw, what is the probability that more than three throws (ie. at least four) will be necessary? ANS: 197/198

Suppose that an ace turns up on an even throw. What is the probability that it turned up on the second throw?

(If you think this is too hard mathematics, do it as a computer assignment!)

[ 5 A deck of 52 cards has 4 suits and 13 values per suit.

1. Write a program simulating the hand of 5 cards.
2. Use your program to answer the following questions:

   1. How often does a pair occur?
   2. How often does a two-pair occur?
   3. How often does a three of a kind (three of same value and two different) occur?
   5. How often does a four of a kind occur?
   6. How often does a full house (2+3) occur?
   7. How often does a straight (five cards in a row, not all same suit) occur?
   8. How often does a flush (five cards of one suit, not in order) occur?
   9. How often does a straight flush (five cards in a row all same suit) occur?

   If you think this is too difficult on a computer, compute the probabilities by hand.

[ 10 A math teacher in a certain school likes to give multiple choice tests, grade them as either right, or wrong, and then lets the students to go over the test and correct the ones they got wrong. This gives them two chances to get a problem right, and the chance of getting a question right increases even if the student just guesses the answers. Suppose a student simply guessed the first time, got the corrections, and guessed differently the second time on the wrong answers. How much his grade improves?

[ 6 This is the expanded version of Exercise 1.5.1. In a group of n people, how often at least two have the same birthday?

1. Find the formula for the probability p(n) assuming 365 days per year, and equally likely birthdays.
2. Compute the probabilities for n = 20, 30, 40, 50
3. Write a simulation program, and verify if the simulation agrees with the theoretical answers.
4. Modify the simulation program to allow for not equal birthdays. Assume January, February, March are less likely then the other days of the year. Change the parameters, and verify how the probabilities p(n) change as you depart from the uniform probabilities. If the change of p(n) is of the magnitude comparable to the simulation accuracy, clearly it is irrelevant.
5. A randomly selected person has chance 1/4 to be born on a leap year. How does this affect the answers?

Chapter 2
Random variables (continued)

Tree species are not distributed at random but are associated with special habitats.
The Auborn Society Field Guide to North American Trees.

Intuitively, random variables are numerical quantities measured in an experiment. The concept¹¹ is the core of probability theory; it leads outside of elementary probability and it touches advanced concepts of integration, function transforms and weak limits.

For convenience random variables are split into three groups: continuous, discrete, and the rest.

2.1  Discrete r. v.

The definition says that X is a discrete r. v. if there is a finite, or countable set \cal V of numbers (values) of X such that p_v = Pr(X = v) > 0 and å_{v Î \cal V} p_v = 1. The function f(v) = p_v is then called the probability mass function of X. For completeness, the domain of the probability mass function is often extended to all x Î IR (or to x Î IR^d in the multivariate case) by f(x) = Pr(X = x).

It is easy to see that if f is a function which satisfies two natural conditions:

For discrete r. v. the cumulative distribution function (12) plays lesser role. It is a discontinuous function given by the expression

2.1.1  Examples of discrete r. v.

Name	Values	Probabilities	Symbol	Parameters
Binomial	0,...,n	Pr(X = k) = (nk)pk(1-p)n-k	Bin(n,p)	0 £ p £ 1, n Î IN
Poisson	ZZ\scriptscriptstyle +	Pr(X = k) = e-l[(lk)/ k!]	Poiss(l)	l > 0
Geometric	ZZ\scriptscriptstyle +	Pr(X = k) = p(1-p)k-1		0 £ p £ 1
Uniform	{x1,...,xk}	Pr(X = xj) = 1/k		k Î IN, x1,...,xk Î IR
Hypergeometric
Negative Binomial

Example 14 Let the random variable X denote the number of heads in three tosses of a fair coin.

Example 15 Let the random variable X denote the score of a randomly selected student on the final exam.

Problem 8 Let N be Poiss(l), and assume N balls are placed at random into n boxes. Find the probability that exactly m of the boxes are empty. ANS: (_mⁿ)e^-lm/n(1-e^l/n)^n-m.

2.1.2  Simulations of discrete r. v.

Other methods are also available for the distributions from Table 2.1. For example, program TOSSCOIN.BAS on page pageref simulates binomial distribution Bin(n=100, p=1/2).

The following exercise provides tools to run more involved simulations. [ 11 Write functions SimOneBin(n,p) and SimOneGeom(p) that will simulate a single occurrence of the Bin(n,p) r. v. and the geometric r. v. The sample usage: PRINT SimOneBin(15,.5) should simulate the number of heads in tossing 15 fair coins.

Also write function SimGeneric(p()) which simulates generic r. v. with values {0,1,...,n} and prescribed probabilities p_k = p(k).

2.2  Continuous r. v.

Since probability satisfies continuity axiom (1), Pr(X Î (a,a+h))® 0 as h® 0 for all a. The main interest in continuous case is that the rate of convergence to 0 is also known, Pr(X Î (a,a+h)) » f(a) h +o(h). Function f(x) in this expansion is called the density function.

In terms of the cumulative distribution function 12), the probability is Pr(X Î (a,a+h)) = F(a+h)-F(a), and thus f(a) = lim_{h® 0}[(F(a+h)-F(a))/ h] = F¢(a) is the derivative of the cumulative distribution function F. Therefore when the Fundamental Theorem of Calculus can be invoked (say, when f is piecewise continuous)

Definition 5 Random variable X is (absolutely) continuous, if there is a function f such that Pr(X Î U) = ò_U f(x)dx. Function f is called the probability density function of X.

It is known that if f is a function which satisfies two natural conditions:

It is convenient to use the heuristic probability density function in continuous case corresponds to probability mass function in discrete case, and that expressions that involve in discrete case sums are replaced by integrals, compare (17) and (16).

2.2.1  Examples of continuous r. v.

Name	Range	Density	Symbol	Parameters
Normal	-¥ < x < ¥	f(x) = [1/(Ö{2p}s)]exp-[((x-m)2)/(2s2)]	N(m,s)	m Î IR,s > 0
Exponential	x > 0	f(x) = le-lx		l > 0
Uniform	a < x < b	f(x) = [1/( b-a)]	U(a,b)	a < b real
Gamma	x > 0	f(x) = Cxa-1e-x/b	Gamma(a,b)	a > 0,b > 0,C = [1/( baG(a))]
Weibull

Example 16 A dart is thrown at a circular dart board of radius 6. Let X denote the distance of the dart from the center. Assuming uniform assignment of probability (3), the density of X is

f(x) = {

[x/ 18]	if 0 £ x £ 6
0	otherwise

Problem 9 Referring to Exercise 1.2.2, let X,Y denote the arrival times of the two drivers at the intersection. Find the density of the time lapse |X-Y| between their arrivals.

2.2.2  Histograms

To create a useful histogram, split the range into the finite number of intervals. Then graph over each interval the rectangle with the area equal to the observed frequency. The number, and positioning of intervals depends on the amount of data available, and personal preference.

2.2.3  Simulations of continuous r. v.

The generic method for simulating a continuous random variable is similar to the method used in the discrete case. Namely, we take X = f(U) with the suitable function f.

To find f assume it is increasing and thus invertible with inverse g. Then Pr(f(U) < x) = Pr(U < g(x)) = g(x). This completes the generic prescription: take as f(x) the inverse¹² of the cumulative distribution functions F(x) = Pr(X £ x).

This method of simulation is quite effective if the inverse of F can be found analytically. It becomes slow when the inverse (or, worse still, cumulative distribution function F itself) is computed by numerical procedures. Since this is the case of the normal distribution, special methods are used to simulate the normal distribution.

Example 17 To simulate X which is exponential with parameter l, use X = -[1/( l)]lnU.

2.3  Expected values

Expected value captures the intuition of the average of a random quantity. It is also this intuition that leads to estimating the expected value by averaging the outcomes of simulations. Example 18 If X has values x₁,...,x_n with equal probability, then EX = [`x] is the arithmetic mean of x₁,...,x_n.

Simulationsare oten used to get answers that are too difficult to find analytically. The following exercise can be answered by simulation if you figure out how to shuffle cards from a deck at random (Exercise 1.4). [ 12 What is the expected number of cards which must be turned over in order to see each of one suit.

Example 19 If X takes two values a < b and Pr (X = a) = p, then EX = pa+(1-p)b is the number in the closed interval [a,b].

Definition 7 For continuous random variable X the expected value EX is given by EX = ò_IRxf(x) dx, provided the integral converges.

Name	Probability distribution	EX
Normal	f(x) = [1/(Ö{2p}s)]exp-[((x-m)2)/( 2s2)]	m
Exponential	f(x) = le-lx	[1/( l)] > 0
Uniform	f(x) = [1/( b-a)]	1/2(a+b) real
Gamma
Weibull
Binomial	Pr(X = k) = (nk)pk(1-p)n-k	np
Poisson	Pr(X = k) = e-l[(lk)/ k!]	l
Geometric	Pr(X = k) = p(1-p)k-1	1/p
Hypergeometric
Negative Binomial

Problem 10 Compute EX for the entries in Table 2.3.

[ 13 Simulate EX for the entries in Table 2.3 (except normal) for different values of parameters involved.

Problem 11 [Exercise] Referring to Example 2.2.1, what is the average distance of a dart from the center?ANS: 4.

If you chose to do the simulations, you should perhaps notice that it is rather difficult to decide how many simulations to take in order to achieve the desired accuracy. Typically, you need to increase the size of a simulation four times to half the error.

Another point to keep in mind is that simulations do return numbers. But from numbers alone it is difficlut to see how parameters of the model change it, and this is a more interesting question.

2.3.1  Tail integration formulas

Proof. To simplify the proof and expose the main idea more clearly, consider the discrete case with a finite number of values. Similarly, to avoid technicalities we consider continuous case with bounded range only.

Discrete case:

EX = x₁p₁+x₂p₂+...+x_np_n = x₁(p₁+p₂+...+p_n-( p₂+...+p_n))+x₂( p₂+...+p_n-( p₃+...+p_n))+...+x_n-1( p_n-1+p_n -p_n)+x_np_n = x₁( p₁+p₂+...+p_n)+(x₂-x₁)( p₂+...+p_n-(p₃+...+p_n))+...+( x_n-1-x_n-2)( p_n-1+p_n)+(x_n-x_n-1)p_n The latter is ò₀^x_n Pr(X > t)  dt.

Continuous case: Let f denote the density and F be the cumulative distribution function. Integrating by parts EX = ò₀^b x f(x)  dx = b-ò₀^bF(x) dx = ò₀^b(1-F(x) dx. ^[¯]
If X is discrete integer valued X Î {0,1,...}, then (21) can be written as

It is natural to define EX for more general r. v. in the same way through formula (21). Write X = X⁺-X^- to decompose X into its non-negative, and non-positive parts, and then define EX = ò₀^¥ P(X > t) dt-ò₀^¥ P(X < -t) dt. Clearly, if one of the integrals diverges, EX is not defined.

Example 20 Suppose X is exponential with the density f(x) = e^-x. Let Y be X truncated at level 3. That is

Y = {

X	if X £ 3
3	if X ³ 3

Clearly, Y is not continuous, as Pr(Y = 3) = Pr(X > 3) = e^-3 > 0. On the other hand, Y is not discrete as it takes uncountable number of values; in fact all the numbers between 0 and 3 are possible. The definitions of the expected value we gave do not apply, but (21) can be used to show that E(Y) = 1-e^-3.

Indeed, Pr(Y > t) = Pr(X > t) = e^-t for 0 < t < 3, and Pr(Y > t) = 0 for t > 3.

Example 21 To determine the mean of the geometric distribution we can either compute the sum på_{n = 1}^¥ n (1-p)^n-1, or use (21) and find the easier sum å_{n = 0}^¥ (1-p)ⁿ.

2.3.2  Chebyshev-Markov inequality

[ 1 If U ³ 0 then

Indeed, by (21) we have EU = ò₀^¥ Pr(U > x) dx ³ ò₀^tPr(U > x) dx ³ tP(|X| > t).

Problem 12 Suppose U is uniform U(0,1). Then Pr(U > t) £ [1/ 2t]. This means that 1-t < [1/ 2t].

1. Is the above inequality ``sharp"? (Graph both curves).
2. Is (21) sharp? That is, given t > 0, is there an X > 0 such that equality occurs?

2.4  Expected values and simulations

Expected value EX is approximated without much difficulty¹³ on a computer by averaging large number of independent trials.

Example 22 To find the average of the uniform U(0,1) distribution, take ¹/_n(U₁+...+U_n).

This is the basis of Monte-Carlo methods , which is a family of related probabilistic methods for computing the integrals. To find ò_a^b f(x) dx we simulate X_j = f(a+bU_j). The average ¹/_nå_{j = 1}ⁿX_j approximates [1/( b-a)]ò_a^b f(u) du for large n.

The variance of the n-th approximation is of the order n^-1/2, which is worse than the trapezoidal rule for smooth functions. In return the approximation is insensitive to the smoothness of the integrands, and also to the dimension of the integral. Monte Carlo methods can be used effectively for multiple integrals over irregular domains.

[ 14 Use Monte Carlo method to approximate p = ò_-1¹2[Ö(1-x²)]  dx. (You may want to compare the output with numerical procedures described in Section .)

Another method of similar nature is to pick points (X,Y) at random from the rectangle containing the graph of f and check if Y < f(X) holds. The proportion of ``successes" approximates the proportion of the area under the graph of f.

[ 15 Approximate p/4 by selecting points (X,Y) at random from the unit square, and checking if X²+Y² < 1.

The following sample program computes numerically double integral òò_Ucos(10x+20y) dxdy over a circle x²+y² = 1. The only conceptual difficulty as compared to single integrals is how to select points at random from the unit disk. This is done by picking points from a bigger square and discarding those that didn't make it. Can you do this integral analytically? Or by another numerical procedure?

PROGRAM dblint.bas
'

'declarations
DECLARE FUNCTION Integrand! (X!, Y!)
DECLARE FUNCTION InDomain! (X!, Y!)

CONST True = -1

' simulation loop
NumTrials = 10000
FOR j = 1 TO NumTrials
        'select random points from the square [-1,1]x[-1,1]
        X = 2 * RND(1) - 1
        Y = 2 * RND(1) - 1
        'check if this is in the domain
        IF InDomain(X, Y) THEN
                NumTested = NumTested + 1
                Sum = Sum + Integrand(X, Y)
                Var = Var + Integrand(X, Y) ^ 2
        END IF
NEXT j


'Print the answer
PRINT "Examined "; NumTested; " random points"

IF NumTested = 0 THEN END 'nothing found
N = NumTested
PRINT "The integral is approximately "; Sum / N
PRINT "With 95% confidence  the error is less than "; 1.96 * SQR(Var / N - (Sum / N) ^ 2) / SQR(N)

END

FUNCTION InDomain (X, Y)
'This function checks if $x,y$ is in the integration domain
'The definition of the domain can be easily modified here, including
'more complicated domains
IF X ^ 2 + Y ^ 2 < 1 THEN InDomain = True
END FUNCTION

FUNCTION Integrand (X, Y)
'This is the function to be integrated

Integrand = COS(10 * X + 20 * Y)
'

It is important to have some idea about how accurate the answer is. The program uses the to estimate the magnitude of the error. A less sharp (and thus safer to use!) error estimate can be obtained from (23)

[ 7 There are two natural choices to estimate by simulation events of small probability. We can pick a large sample size n, run the simulation experiment and hopefully get several data points. The outcome X of such an experiment is a binomial random variable with unknown probability p of success, and we would estimate p » X/n. The trouble is that if we don't know how small the chances are, we might get none, estimating probability to be zero. Or we could run the experiment until we get the prescribed number of successes. The observation would then consist of a set of geometric random variables T₁,..., T_k with unknown mean ET = 1/p. We could then estimate p = 1/ET by taking the inverse of the arithmetic mean of T_j. In this approach we are guaranteed to get some observations, as long as p ¹ 0.

The question is Which of the methods would you recommend to use? (Of course, you would recommend a better method, but what the word ``better" might mean here?)

2.5  Joint distributions

Often an experiment involves measuring two or more random numbers, say X and Y. The fact that we know the distribution of X, and the distribution of Y separately doesn't determine probabilities of events that involve both X and Y simultaneously. Example 23 Suppose

It is clear that if X is discrete, and Y is discrete, then (X,Y) is an IR² valued discrete r. v. That is, the values of the pair are countable. Probabilities Pr(X = x, Y = y) are called the joint distribution. Corresponding Pr(X = x) and Pr(Y = y) are the so called marginals. Example 2.5 points out that marginals do not determine joint probabilities uniquely. But if we know the joint probabilities then we can compute the marginals, eg Pr(X = x) = å_yPr(X = x, Y = y).

In contrast to the discrete case, joint continuity can not be recognized from the continuity of the components, and requires full definition.

Definition 8 Let X = (X₁,...,X_n). Random variables X₁,...,X_n are jointly (absolutely) continuous, if there is a function f such that

Example 24 Suppose X,Y have uniform distribution in the unit disk. Then the joint density is f(x,y) = 1/p for x²+y² £ 1 and the density of X is f_X(x) = [2/( p)][Ö(1-x²)].

The relation between probabilities and the density is

2.5.1  Independent random variables

(Similarly we define the stochastic independence of random vectors).

We say that X₁,...,X_n are independent identically distributed (i. i. d) if (26) holds and Pr(X_i Î U) = Pr(X_j Î U) for all Borel U Ì IR. [ 2 If X,Y are discrete with the probability mass function f(x,y), then independence of X,Y is equivalent to f(x,y) = f_X(x)f_Y(y).

If X,Y are continuous with the joint density f(x,y) then independence of X,Y is equivalent to f(x,y) = f_X(x)f_Y(y).

Independence is often part of the model. Independence allows to determine joint distributions from marginals. Thus each independent random variable can be analyzed separately, and then more complex questions can be answered. From the mathematical perspective, under independence we can determine joint distributions if we know marginals. Example 25 Suppose X is binomial Bin(n,p) and Y is Poiss(l). If X,Y are independent, then the joint probability mass function of X,Y is given by f(x,y) = (_xⁿ)pⁿ(1-p)^n-xe^-ll^y/y! (or 0).

Example 26 [Example 1.2.2 continued] Two drivers arrive at an intersection between 8:00 and 8:01 every day. Their arrival times are independent random variables. Indeed, using formula (25) and elementary area computation, Pr(X < a,Y < b) = ab for 0 < a,b < 1.

2.6  Functions of r. v.

Sums and linear combinations, medians, maxima, and minima are perhaps the most often encountered functions of multidimensional random variables. Methods to compute the distribution, or the expected value of such a function are of considerable practical significance.

2.7  Moments of functions

If X,Y are discrete and EZ exists, then

If X,Y are jointly continuous then Z might be continuous, discrete, or say of mixed type. Regardless of the case

In particular the expected value is linear

The fact that expected value is linear provides a simple method of computing some otherwise difficult sums. Example 27 Suppose X is Binomial Bin(n,p). Then X = X₁+...+X_n, where X_j is the number of successes in j-th trial. Clearly each X_j is 0, or 1, and EX_j = p.

[ 16 [Example 1.2.2 continued] Two drivers arrive at an intersection between 8:00 and 8:01 every day. On average how much time lapses between their arrivals?

Definition 10 Let m = EX. The variance Var(X) is defined as Var(X) = E(X-m)².

Notice that

The standard deviation is s = [ÖVar(X)].

Name	Probability distribution	Var(X)
Normal	f(x) = [1/(Ö{2p}s)]exp-[((x-m)2)/( 2s2)]	s2
Exponential	f(x) = le-lx	[1/( l2)]
Uniform	f(x) = [1/( b-a)]	[1/ 12](b-a)2
Gamma
Weibull
Binomial	Pr(X = k) = (nk)pk(1-p)n-k	np(1-p)
Poisson	Pr(X = k) = e-l[(lk)/ k!]	l
Geometric	Pr(X = k) = p(1-p)k-1	[1/( p2)]
Hypergeometric
Negative Binomial

Problem 13 Compute variances Var(X) for the entries in Table 2.4. (Some of these are a real challenge to your computational skills, so you may safely give up. Another method will make it easier in Chapter ).

Since Var(X) = EX²-(EX)² ³ 0, the following inequality follows

Definition 11 The covariance of random variables X,Y with expected values m_X, m_Y is defined as cov(X,Y) = E(X-m_X)(Y-m_Y).

Clearly Cov(X,X) = Var(X) and Var(X+Y) = Var(X)+Var(Y)+2cov(X,Y).

Theorem 3 If X,Y are independent, then Var(X+Y) = Var(X)+Var(Y).

Example 28 let X₁,X₂,...,X_n be independent {0,1}-valued random variables, and suppose that Pr(X_j = 1) = p. Then Var(å_{j = 1}ⁿX_j) = np(1-p). What is the distribution of å_{j = 1}ⁿX_j?

Tail integration formula revisited

Example 29 If X > 0 then Ee^X = 1+ò₀^¥ e^t P(X > t) dt

Problem 14 Show that EX² = ò₀^¥ t P(|X| > t) dt.

Generalize this formula to E|X|^p.

Chebyshev-Markov inequality

Special cases of Chebyshev's inequality (23) are:

Chebyshev's inequality is one reason we often strive for small average quadratic error. If E|X-X₀|² < e then we can be sure that Pr(|X-X₀| > ³Ö{e}) < ³Ö{e}. The following may be used (with some caution) in computer programs to asses error in estimating probabilities by sampling.

Example 30 If X is Bin(n,p) then Pr(|^X/_n-p| > [1/( ³Ön)]) £ [1/( 4³Ön)]

[ 17 Run a simulation of the event that you know probability of, printing out the error estimate.

Do the same for the event that you don't know the probability analytically. (Use programs written for previous exercises)

2.7.1  Simulations

From now on, in the output of your simulation programs you should provide some error estimates.

2.8  Application: scheduling

2.8.1  Deterministic scheduling problem

As an example of simple scheduling problem consider the following. Example 31 Suppose that we want to bake a batch of chocolate-chip cookies. The tasks and their (estimated) times are:

• [T1] Bake at 350F (40 min)
• [T2] Make batter (5 min)
• [T3] Pre-heat oven to 350F (10 min)
• [T4] Find and grease pan (2 min)
• [T5] Find a cookie-tray to serve cookies (2 min)
• [T6] Take cookies out, cool and serve (5 min)

The dependency graph is quite obvious here; for instance, we cannot start baking before batter is ready. What is the shortest time we can eat the cookies?

2.8.2  Stochastic scheduling problem

[ 18 Suppose for the sake of this exercise that the numbers presented in the cookie-baking example are just the average values of the exponential random variables.

• What will be the average completion time then? Will it be larger, smaller, or about equal to the previous answer?
• How often will we finish the process before the previously estimated (deterministic) time?

2.8.3  More scheduling questions

2.9  Distributions of functions

Sums of discrete r. v.

Suppose X,Y are discrete and f(x,y) = Pr(X = x, Y = y) is their joint probability mass function. Then Z = X+Y is discrete with values z = x+y. Therefore

For independent random variables X, Y this takes a slightly simpler form.

Formula (36) can be used to prove the so called summation formulas. Theorem 4 If X, Y are independent Binomial with the same parameter p, ie. X is Bin(n,p) and Y is Bin(m, p), then X+Y is Binomial Bin(n+m,p).

If X, Y are independent Poisson Poiss(l_X) and Poiss(l_Y), then X+Y is Poisson Poiss(l_X+l_Y).

Sums of continuous r. v.

One can prove that if X,Y are independent and continuous than X+Y is continuous with the density

Formula (37) defines the convolution f_X*f_Y. It can be used to prove the so called summation formulas. Theorem 5 If X, Y are independent Normal, then X+Y is Normal.

If X, Y are independent Gamma with the same parameter b, then X+Y is Gamma(a_X+a_Y,b).

Example 33 [Example 1.2.2 continued] Two drivers arrive at an intersection between 8:00 and 8:01 every day. What is the density of the time that lapsed between their arrivals?

Example 34 Suppose X,Y are independent U(0,1). The density of Z = X+Y is

f(z) = {

z	if 0 £ z £ 1
2-z	if 1 £ z £ 2

.

Minima, maxima

Suppose X,Y are independent. If U = min{X,Y} then Pr(U > t) = Pr(X > t)Pr(Y > t). Therefore the reliability function of U can be computed from the two given ones.

Example 35 If X, Y are independent exponential, then U = min{X,Y} is exponential.

If U = max{X,Y} then Pr(U < t) = Pr(X < t)Pr(Y < t). Therefore the cumulative distribution function of U can be computed from the two given ones.

Example 36 Suppose X,Y are independent uniform U(0,1). Then U = max{X,Y} has the density f(u) = 2u for 0 < u < 1.

Problem 15 Let U₁,..., U_n be independent uniform U(0,1). Find the density of

• X = min_j U_j
• Y = max_j U_j

Problem 16 If X,Y are independent exponential random variables with parameters l, m, show that Pr(X < Y) = [(l)/( l+m)].

Order statistics

Let R₁,..., R_n be the corresponding order statistics. This means that at the end of each experiment we re-arrange the numbers X₁,... X_n into the increasing sequence R₁,..., R_n. This means that R₁ = min_j X_j is the smallest, R₂ = max_imin_{j ¹ i} X_j is the second largest, etc.

The density of R_k can be found by the following method. In order for the inequality R_k > x to hold, there must be at least k values among the X_j above level x. Since X_j are independent and have the same probability p = Pr(X_j > x) of ``success" in crossing over the x-level, this means that Pr(R_k > x) is given by the binomial formula with n trials and probability of success p = 1-G(x).

When the derivative is taken, the sum collapses into just one term, giving the elegant answer r_k(x) = [n!/( (k-1)!(n-k)!)](G(x))^k(1-G(x))^n-kg(x).

2.10  L₂-spaces

The L₂ norm is

Notice that ||X-EX||₂ is just another notation for the standard deviation. Thus standard deviation is the L₂ distance of X from a constant.

We say that X_n converges to X in L₂, if ||X_n-X||₂® 0 as n® ¥. We shall also use the phrase sequence X_n converges to X in mean-square. An example of the latter is Theorem .

Several useful inequalities are collected in the following¹⁴. Theorem 6 For all square-integrable X,Y

• Cauchy-Schwarz inequality:

  EXY £ ||X||2||Y||2.

  (39)

• Jensen's inequality:

  E|X| £ E||X||2.

  (40)

• Triangle inequality:

  ||X+Y||2 £ ||X||2+||Y||2.

  (41)

Proof. Proof of (39): Quadratic function f(t) = E|X+tY|² = EX²+2tEXY+t²EY² is non-negative for all t. Therefore its determinant D ³ 0. (Compute D to finish the proof.)

Proof of (40): Use (39) with Y = 1.

Proof of (41): By (39) E|X+Y|² £ ||X||₂²+||Y||₂²+2||X||₂||Y||₂. ^[¯]

2.11  Correlation coefficient

Random variables with r = 0 are called uncorrelated. Correlation coefficient close to one of the extremes ±1 means that there is a strong linear relation between X,Y; this is stated more precisely in ().

Theorem 2.7 states that independent random variables with finite variances are uncorrelated.

Problem 17 Give an example of dependent random variables that are uncorrelated.

2.11.1  Best linear approximation

Suppose we would like to approximate random variable Y by another quantity X that is perhaps better accessible. Of all the possible ways to do it, linear function Y » mX+b is perhaps the simplest. Of all such linear functions, we now want to pick the best. In a single experiment, the error is |Y-mX-b|. We could minimize the average empirical error over many experiments ¹/_nå_j|Y_j-mX_j-b|. This approximates the average error E|Y-mX-b|. Let us agree to measure the error of the approximation by a quadratic error E(Y-mX-b)² instead. (This choice leads to simpler mathematics.)

Question: For what values of m, b the error E(Y-mX-b)² is the smallest? When is it 0?

Let H(m,b) = E(Y-mX-b)². Clearly, H is a quadratic function of m,b Î IR. The unique minimum is determined by the system of equations

The answer is m = [cov(X,Y)/ Var X], b = EY-mEX, and the minimal error is

2.12  Application: length of a random chain

A chain in the x,y-plane consists of n links each of unit length. The angle between two consecutive links is ±a, where a > 0 is a constant. Assume the sign is taken and random, with probability ¹/₂ for each. Let L_n be the distance from the beginning to the end of the chain. The angle between the k-th link and the positive x-axis is a random variable S_k-1, where we may assume (why?) S₀ = 0 and S_k = S_k-1+x_ka, where x = ±1 with probability ¹/₂. The following steps determine the average length of the random chain.

1. L_n² = (å_{k = 0}^n-1cosS_k)²+(å_{k = 0}^n-1sinS_k)².
2. EcosS_n = cosⁿ a
3. EsinS_n = 0
4. EcosS_m cosS_n = cos^n-maEcos² S_m for m < n
5. EsinS_m sinS_n = cos^n-maEsin² S_m for m < n
6. EL_n²-L_n-1² = 1+2cosa[(1-cos^n-1a)/( 1-cosa)]

7. EL_n² = n [(1+cosa)/( 1-cosa)] -2cosa[(1-cosⁿa)/( (1-cosa)²)]

2.13  Conditional expectations

2.13.1  Conditional distributions

Definition 12 The conditional expectation E{X|Y = y} is defined as åxPr(X = x|Y = y) in discrete case, and as ò_IR x f(x|Y = y)  dx in the continuous case. One can show that the expected values exist, when E|X| < ¥.

Example 37 A game consists of tossing a die. If the face value on the die is X then a coin is tossed X times. Let Y be the number of heads. Then E(Y|X = x) = ¹/₂x.

2.13.2  Conditional expectations as random variables

Example 38 Suppose Y is a discrete with different values on the events A₁, A₂, ¼, A_n which form a non-degenerate disjoint partition of the probability space W. Then

Example 39 Suppose that f(x, y) is the joint density with respect to the Lebesgue measure on IR² of the bivariate random variable (X, Y) and let f_Y(y) ¹ 0 be the (marginal) density of Y. Put f(x|y) = f(x, y)/f_Y(y). Then E{X|Y} = h(Y), where h(y) = ò_-¥^¥ x f(x|y) dx.

Total probability formula for conditional expectations is as follows.

2.13.3  Conditional expectations (continued)

Similar question can be solved for N having a Poisson distribution.

Example 42 What is the distribution of a geometric sum of i. i. d. exponential r. v.?

Example 43 Stock market fluctuations can be modelled by Z = x₁+...+x_N, where N, the number of transactions, is Poisson(l) and x are normal N(0,s). There is no explicit formula for the density of Z, but there is one for the moment generating function. Thus Chebyshev inequality gives bounds of the form Pr(Z > t) £ exp....

2.14  Best non-linear approximations

Theorem gives geometric interpretation of the conditional expectation E{·|Z}; for square integrable functions E{.|Z} is just the orthogonal projection of the Banach (normed) space L₂ onto its closed subspace L₂(Z), consisting of all 2-integrable random variables of the form f(Z).

Theorem 7 For square integrable Y the quadratic error E(Y-h(X))² among all square integrable functions h(X) is the smallest if h(x) = E(Y|X = x).

Theorem 2.14 implies that the linear approximation from Section is usually less accurate. In Chapter we shall see that linear approximation is the best one can get in the all important normal case. Even in non-normal case linear approximations offer quick solutions based on simple second order statistics. In contrast, the non-linear approximations require elaborate numerical schemes to process the empirical data.

2.15  Lack of memory

Suppose T represents a failure time of some device. If the device is working at time t, then the probability of surviving additional s seconds is Pr(T > t+s|T > t). For a device that doesn't exhibit aging this probability should be the same as for the brand new device.

Problem 19 Show that geometric T satisfies (47) for integer t,s.

Equation (47) implies Pr(T > t+s) = Pr(T > t)Pr(T > s), an equation that can be solved. Theorem 8 If T > 0 satisfies (47) for all t,s > 0, then T is exponential.

Proof. The tail distribution function N(x) = P(T > x) satisfies equation

Formula (48) implies that for all integer n and all x ³ 0

This shows that for rational q > 0

2.16  Intensity of failures

Let T > 0 be a continuous r. v. interpreted as a failure time of a certain device. If the device is in operational condition at time t, then the probability that it will fail immediately afterwards may be assumed negligible. The probability of failing within h units of time is Pr(T < t+h|T > t). The failure rate at time t is defined as

Example 44 If T is exponential then the failure rate is constant.

A family of failure rates that exhibit interesting aging patterns is provided by the family of power functions l(t) = t^a.

Theorem 9 If T is continuous with failure rate l(t) = t^a, where a > 0 then T has the Weibull density:

2.17  Poisson approximation

Proof. Rewrite the expression (_kⁿ)[(l^k)/( n^k)](1-[(l)/ n])^n-k = l^k/k!(1-l/n)ⁿ/(1-l/n)^kÕ_{j = 0}^k(1-j/n). ^[¯]
Example 45 How many raisins should a cookie have on average so that no more than one cookie in a hundred has no raisins? So that no more than one cookie in a thousand has no raisins?

Example 46 A bag of chocolate chip cookies has 50 cookies. The manufacture claims there are a 1,000 chips in a bag. Is it likely to find a cookie with 15 or less chips in such a bag?

2.18  Questions

From numbers 1, ... n, select at random k > n numbers. On average, how many of these numbers repeat? (and hence should be thrown out)

Problem 21 The performance of the algorithm for selecting a random permutation in GetPermutation SUB of RANDTOUR.BAS can be estimated by analyzing the following ``worst-case" scenario.

When the algorithm attempts to select last of the random numbers 1, ... n, then

1. What is the probability if will find the ``right number" on first attempt?
2. How many attempts on average does the algorithm take to find the last random number?

[ 19 What is the probability that in the group of 45 people one can find two born on the same day of the year? (Compare Example 1.5.1).

Problem 22 For a group of n person, find the expected number of days of the year which are birthdays of exactly k people. (Assume 365 days in a year and that all birthdays are equally likely.)

Problem 23 A large number N of people are subject to a blood test. The test can be administered in one of the two ways:

(i) Each person can be tested separately (N tests are needed).

(ii) The blood sample of k people can be pooled (mixed) and analyzed. If the test is positive, each of the k people must be tested separately and in all k+1 tests are then required for k people.

Assume the probability p that the test is positive is the same for all people and that the test results for different people are stochastically independent.

1. What is the probability that the test for a pooled sample of k people is positive?

2. What is the expected number of tests necessary under plan (ii)?

3. Find the equation for the value of k which will minimize the expected number of tests under plan (ii).
4. Show that the optimal k is close to [1/( Öp)] and hence that the minimum expected number of tests is on average about [2N/( Öp)].

[ 20 Solve Exercise 1.2.2 on page pageref assuming that the arrival times are exponential rather than uniform. Assume independence.

[ 21 There are 3 stop-lights spaced within 1km of each other and operating asynchronously. (They are reset at midnight.) Assuming each is red for 1 minute and then green for one minute, what is the average time to pass through the three lights by a car that can instantaneously accelerate to 60km/h.

This exercise can be developed into a simulation project that may address some of the following questions

• How does the speed change with the number of lights?
• How does the answer change if a car has finite acceleration?
• Are the answers different, if each green light lasts random amount of time, 1min on average?
• How to model/simulate more than one car?
• Can you simulate car traffic on a square grid with stop-lights at intersections?

More theoretical questions

Problem 25 Let X ³ 0 be a random variable and suppose that for every 0 < q < 1 there is T = T(q) such that

Problem 26 If x, U are independent, Pr(x = 0) = Pr(x+1) = ¹/₂ and U is uniform U(0,1). What is the distribution of U+x?

Chapter 3
Moment generating functions

3.1  Generating functions

Properties of a sequence {a_n} are often reflected in properties of the generating function h(z) = å_na_nzⁿ.

3.2  Properties

A moment generating function is non-negative, and convex (concave up). The typical example is the moment generating function of the {0,1}-valued random variable.

A linear transformations of X changes the moment generating function by the following formula.

Important properties of moment generating functions are proved in more theoretical probability courses¹⁵.

Theorem 11 (i) The distribution of X is determined uniquely by its moment generating function M(t).

(ii) If X, Y are independent random variables, then M_X+Y(t) = M_X(t) M_Y(t) for all t Î IR.

(iii) M(0) = 1, M¢(0) = EX, M¢¢(0) = EX²

Name	Distribution	Moment generating function
Normal N(0,1)	f(x) = [1/(Ö{2p})]exp-[(x2)/ 2]	M(t) = et2/2
Exponential	f(x) = le-lx	M(t) = [(l)/( l-t)]
Uniform U(-1,1)	f(x) = 1/2 for -1 £ x £ 1	M(t) = 1/tsinht
Gamma	f(x) = 1/G(a)b-a xa-1exp-x/b	M(t) = (1-bt)-a
Binomial	Pr(X = k) = (nk)pk(1-p)n-k	M(t) = (1-p+pet)n
Poisson	Pr(X = k) = e-l[(lk)/ k!]	M(t) = expl(et-1)
Geometric	Pr(X = k) = p(1-p)k-1	M(t) = [(pet)/( 1-(1-p)et)]

Problem 27 Find moment generating functions for each of the entries in Table 3.1.

Problem 28 Use moment generating functions to compute EX, Var(X) for each of the entries in Table 2.4.

Problem 29 Prove the summation formulas stated in Theorems 36 and 37

For a d-dimensional random variable X = (X₁, ¼, X_d) the moment generating function M_X: IR^d® \sf CC is defined by M_X(t) = Eexp(t·X), where the dot denotes the dot (scalar) product, ie. x·y = åx_ky_k. For a pair of real valued random variables X, Y, we also write M(t, s) = M_{(X, Y)}((t, s)) and we call M(t, s) the joint moment generating function of X and Y.

The following is the multi-dimensional version of Theorem 3.2.

Theorem 12 (i) The distribution of X is determined uniquely by its moment generating function M(t).

(ii) If X, Y are independent IR^d-valued random variables, then

3.2.1  Probability generating functions

3.3  Characteristic functions

The major nuisance in using the moement generating functions is the fact that the moment generating functions may not exist, when the definiting integral diverges.

For this reason it is more preferable to use an expression that is always bounded, and yet has the same convenient algebraic properties. The natural candidate is

For symmetric random variables complex numbers can be avoided at the expense of trigonometric identities. Example 47 If X is -1,1 valued, Pr(X = 1) = ¹/₂, then f(t) = cost.

Example 48 The characteristic function of the normal N(0,1) distribution is f(t) = e^-t2/2.

3.4  Questions

Problem 31 Suppose the probability p_n that a family has exactly n children is apⁿ when n ³ 1 and suppose p₀ = 1-a[p/( 1-p)]. (Notice that this is a constaint on the admissible values of a, p since p₀ ³ 0.

Suppose that all distributions of the sexes for n children are equally likely. Find the probability that a family has exactly k girls.
Hint: The answer is at first as the infinite series. To find its sum, use generating function, or negative binomial expansion: (1+x)^-k = 1+[(k+1)/ 1!]x+[((k+1)k+2)/ 2!]x²+....

ANS:

[(2ap^k)/( (2-p)^k+1)].

Problem 32 Show that if X ³ 0 is a random variable such that

Problem 33 Show that if Eexp(lX²) = C < ¥ for some a > 0, then

Problem 34 Prove that function f(t): = Emax{X, t} determines uniquely the distribution of an integrable random variable X in each of the following cases:

• [(a)] If X is discrete.

• [(b)] If X has continuous density.

Problem 35 Let p > 2. Show that exp(|t|^p) is not a moment generating function.

Chapter 4
Normal distribution

Next to a stream in a forest you see a small tree with tiny, bell-shaped, white flowers in dropping clusters.
The Auborn Society Field Guide to North American Trees.

The predominance of normal distribution is often explained by the . This theorem asserts that under fairly general conditions the distribution of the sum of many independent components is approximately normal. In this chapter we give another reason why normal distribution might occur. The usual definition of the standard normal variable Z specifies its density f(x) = [1/( Ö{2p})]e^{-[(x2)/ 2]}, see Figure 2.1 on page pageref. In general, the so called N(m, s) density is given by

In multivariate case it is more convenient to use moment generating functions directly. For consistency we shall therefore adopt the following definition.

Definition 13 A real valued random variable X has the normal N(m, s) distribution if its moment generating function has the form

From Theorem 3.2 one can check by taking the derivatives that m = EX and s² = Var(X). Using (53) it is easy to see that every univariate normal X can be written as

If X has given mean m = 123 and given variance s² = 456, for what values of a we have Pr(|X| > a) = .8?

[ 22 [Basketball coach's problem] Collect data on the heights of 25 to 50 randomly selected males. Make a histogram of the data, compute the empirical mean and standard deviation.

• Does it appear that the normal distribution is a good probability distribution for these heights?
• There are about 200,000 males in Cincinnati area. Assuming normal distribution of heights with the mean and variance as you obtained from the data, estimate the number of males taller than 7¢.

4.1  Herschel's law of errors

``Suppose a ball is dropped from a given height, with the intention that it shall fall on a given mark. Fall as it may, its deviation from the mark is error, and the probability of that error is the unknown function of its square, ie. of the sum of the squares of its deviations in any two rectangular directions. Now, the probability of any deviation depending solely on its magnitude, and not on its direction, it follows that the probability of each of these rectangular deviations must be the same function of its square. And since the observed oblique deviation is equivalent to the two rectangular ones, supposed concurrent, and which are essentially independent of one another, and is, therefore, a compound event of which they are the simple independent constituents, therefore its probability will be the product of their separate probabilities. Thus the form of our unknown function comes to be determined from this condition...''

Ten years after Herschel, the reasoning was repeated by J. C. Maxwell¹⁸. The fact that velocities are normally distributed is sometimes called Maxwell's theorem.

The beauty of the reasoning lies in the fact that the interplay of two very natural assumptions: of independence and of rotation invariance, gives rise to the normal law of errors - the most important distribution in statistics.

Theorem 13 Suppose random variables X, Y have joint probability distribution m(dx, dy) such that

(i) m(·) is invariant under the rotations of IR²;

(ii) X, Y are independent.

Then X, Y are normal.

The following technical lemma asserts that moment generating function exists. [ 1 If X,Y are independent and X+Y, X-Y are independent, then EexpaX < ¥ for all a.

Proof. Consider a real function N(x): = P(|X| ³ x). We shall show that there is x₀ such that

Let X₁, X₂ be the independent copies of X. Inequality (56) follows from the fact that event {|X₁| ³ 2x} implies that either the event {|X₁| ³ 2x}Ç{|X₂| ³ 2x₀}, or the event {|X₁+X₂| ³ 2(x - x₀)}Ç{|X₁ - X₂ | ³ 2(x - x₀)} occurs.

Indeed, let x₀ be such that P(|X₂| ³ 2x₀) £ ¹/₂. If |X₁| ³ 2x and |X₂| < 2x₀ then |X₁±X₂| ³ |X₁| - |X₂| ³ 2(x -x₀). Therefore using independence and the trivial bound P(|X₁+X₂| ³ 2a) £ P(|X₁| ³ a)+P(|X₂| ³ a), we obtain

Therefore (take u = v) the second derivative Q¢¢(u) = Q¢¢(0) = const ³ 0. This means M(u) = exp( au+bu²). ^[¯]

4.2  Bivariate Normal distribution

If EX = EY = 0, the moment generating function M(t,s) = Ee^tX+sY is given by M(t,s) = e^{1/₂(s₁2 t2+2tsrs₁s₂+s2t2)}

Clearly s₁² = Var(X), s₂² = Var(Y), r = Cov(X,Y).

When r ¹ ±1 the joint density of X,Y exists and is given by

Example 49 If X,Y are jointly normal with correlation coefficient r ¹ ±1 then the conditional distribution of Y given X is normal.

4.3  Questions

Problem 37 If X,Y are independent normal N(0,1), find the density of X²+Y². (Hint: compute cumulative distribution function, integrating in polar coordinates.)

Problem 38 For jointly normal X,Y show that E(Y|X) = aX+b is linear.

Problem 39 If X,Y are jointly normal then Y-rXs_Y/s_X and X are independent.

Problem 40 If X, Y are jointly normal with variances s_X²,s_Y² and the correlation coefficient r, then X = s_X(g₁ cosq+g₂sinq, Y = s_Y(g₁ sinq+g₂cosq, where g_j are independent N(0,1) and sin2q = r.

Chapter 5
Limit theorems

This is a short chapter on asymptotic behavior of sums and averages of independent observations. Theorem justifies simulations as the means for computing probabilities and expected values. Theorem provides error estimates.

5.1  Stochastic Analysis

Definition 15 X_n® X in probability, if Pr(|X_n-X| > e)® 0 for all e > 0

Example 50 Let X_n be N(0,s = ¹/_n). Then X_n® o in probability.

Definition 16 X_n® X almost surely, if Pr(X_n® X) = 1.

Example 51 Let U be the uniform U(0,1) r. v. Then ¹/_n U® 0 almost surely.

Example 52 Let U be the uniform U(0,1) r. v. et

X_n = {

0	if U > 1/n
1	otherwise

. Then X_n® 0 almost surely.

Definition 17 X_n® X in L₂ (mean square), if E|X_n-X|² ® 0 as n® ¥.

Remark 2 If X_n® X in L₂, then by Chebyshev's inequality X_n® X in probability.

If X_n® X almost surely, then X_n® X in probability.

5.2  Law of large numbers

Theorem 14 Suppose X_k are such that EX_k = m, Var(X_k) = s², and cov(X_i,X_j) = 0 for all i ¹ j. Then ¹/_nå_{j = 1}ⁿ X_j® m in L₂

Proof. To show that ¹/_nå_{j = 1}ⁿ X_j® m in mean square, compute the variance. ^[¯]
[ 1 If X_n is Binomial Bin(n,p) then ¹/_n X_n ® p in probability.

5.2.1  Strong law of large numbers

The following method can be used to prove strong law of large numbers for Binomial r. v.

1. If X is Bin(n,p), use moment generating function to show that E(X-np)⁴ £ Cn²
2. Use Chebyshev's inequality and fourth moments to show that

   å_nPr(|X_n/n-p| > e) < ¥.

3. Use the convergence of the series to show that lim_N®¥Pr(È_{n ³ N} {|X_n/n-p| > e}) = 0.

4. Use the continuity of the probability measure to show that Pr(Ç_NÈ_{n ³ N} {|X_n/n-p| > e}) = 0.

5. Show that with probability one for every rational e > 0 there is N = N(e) such that for all n > N the inequality |X_n/n-p| < e holds. Hint: if Pr(A_e) = 1 for rational e, then Pr(Ç_eA_e) = 1.

5.3  Convergence of distributions

In addition to the types of convergence introduced in Section , we also have the convergence in distribution.

Definition 18 X_n converges to X in distribution, if Ef(X_n)® Ef(X) for all bounded continuous functions f.

Theorem 15 If X_n ® X in distribution, then Pr(X_n Î (a,b))® Pr(X Î (a,b)) for all a < b such that Pr(X = a) = Pr(X = b) = 0.

Theorem 16 If M_Xn(u)® M_X(u) for all u, then X_n® X in distribution.

Theorem 2.17 states convergence in distribution to Poisson limit. Here is a proof that uses moment generating functions.

Proof of Theorem 2.17. M_Xn(u) = (1+p_n(e^u-1))ⁿ® e^l(eu-1). ^[¯]

5.4  Central Limit Theorem

We state first normal approximation to binomial. Theorem 17 If X_n is Binomial Bin(n,p) and 0 < p < 1 is constant, then [(X_n-np)/( [Önpq])]® Z in distribution to N(0,1) random variable Z.

Proof. For p = 1/2 only. The moment generating function of [(X_n-np)/( [Önpq])] = [(2X_n-n)/( Ön)] is

M_n(u) = e^-Önu(¹/₂+¹/₂e^-2u/Ön)ⁿ = ([(e^-u/Ön+e^u/Ön)/ 2])ⁿ = coshⁿ (u/Ön)® e^u2/2. (Can you justify the limit?) ^[¯]

Picture Omitted

Limit theorems are illustrated by the program LIMTHS.BAS. This program graphs empirical proportions [^p]_t as the (random) function of t £ n, makes a histogram, and compares it with the Normal and Poisson histograms. By trying various lengths of path n, one can see the almost sure Law of Large Numbers, and for moderate n see the normal approximation.

Problem 41 The graph of [^p]_n as the function of n as given by LIMTHS.BAS suggests a pair of ``curves" between which the averages are ``squeezed". What are the equations of these curves?

[ 23 Suppose a poll of size n is to be taken, and the actual proportion of the voters supporting an issue in question is p = ¹/₂. Determine the size n such that the observed proportion [^p] = ¹/_nX satisfies Pr([^p] > .8) £ .01.

[ 24 Plot the histogram for a 100 independent Binomial Bin(n = 100,p = .5) random variables.

Theorem 18 If X_j are i.i.d. with EX = m, Var(X) = s² then [(å_{j = 1}ⁿ X_j- nm)/( sÖn)]® Z

Proof. For m = 0, s = 1 only.

The moment generating function is

M_n(u) = (M₁(u/Ön))ⁿ = (1+[u/( Ön)]M¢(0)+[(u²)/ 2n]M¢¢(0)+O([1/( n²)]))ⁿ® e^{u2/2 [¯]}
A lesser known aspect of the central limit theorem is that one can actually simulate paths of certain continuous time processes by taking X_n(t) = [1/( Ön)]å_{k £ nt} x_k, where x_k are independent mean-zero r. v.

The following program uses this to simulate random curves. The program requires a graphics card on the PC. PROGRAM firewalk.bas
'

'This program simulates random walk paths (with uniform incerments)
'reflected at the boundaries of a region
'
'declarations of subs
DECLARE SUB CenterPrint (Text$)

' minimal error handling - graphics card is required
ON ERROR GOTO ErrTrap

CLS
'request good graphics (some cards support SCREEN 7, etc)
SCREEN 9

LOCATE 1, 1 'title
CenterPrint "Path of reflected random walk"
LOCATE 9, 1  'timer location
PRINT "Timer"

scale = 10 '
WINDOW (0, 0)-(scale, scale): VIEW (150, 100)-(600, 300)
LINE (0, 0)-(scale, scale), 10, B': LINE (scale, scale)-(2 * scale, 0), 11, B
FOR j = -4 TO 4
LINE (0, scale / 2 + j)-(scale / 50, scale / 2 + j), 2
NEXT j
X = scale / 2
Y = scale / 2
dead = 0
T = 0
col = 14 * RND(1)
speed = scale / 100
WHILE INKEY$ = "" 'infinite loop until a key is pressed
T = T + 1
X0 = X
Y0 = Y
X = X + (RND(1) - 1 / 2) * speed
Y = Y + (RND(1) - 1 / 2) * speed
IF X < 0 THEN X = -X: col = 14 * RND(1)
IF Y < 0 THEN Y = -Y: col = 14 * RND(1)
IF X > scale THEN X = 2 * scale - X: col = 14 * RND(1)
IF Y > scale THEN Y = 2 * scale - Y: col = 14 * RND(1)
IF X > 2 * scale THEN X = 4 * scale - X: col = 14 * RND(1)
LINE (X0, Y0)-(X, Y), col
LOCATE 10, 1
PRINT T
WEND

END

ErrTrap: 'if there are errors then quit
CLS
PRINT "This error requires graphics card VGA"
PRINT "Error running program"
PRINT "Press any key ...'"
WHILE INKEY$ = ""
WEND
END

SUB CenterPrint (Text$)
' Print text centered in 80 column screen
offset = 41 - LEN(Text$) \ 2
IF offset < 1 THEN offset = 1
LOCATE , offset
PRINT Text$
'

5.5  Limit theorems and simulations

One role of limit theorems is to justify the simulations and provide error estimates. The simulation presented on page pageref uses the Central Limit Theorem to print out the so called 95% confidence interval.

Central limit theorem is also a basis for a fast simulation of the normal distribution, see Section .

5.6  Large deviation bounds

We begin with the bound for binomial distribution.

Recall that the relative entropy is H(q|p) = -qlnq/p - (1-q) ln[(1-q)/( 1-p)].

Theorem 19 Suppose X_n is Binomial Bin(n,p). Then for p¢ ³ p

Proof. Use Chebyshev's inequality (34) and the moment generating function. Pr(X_n ³ n p¢) = Pr(uX_n ³ un p¢) £ e^-np¢u E^uX_n =

( e^-p¢u (1-p+pe^u))ⁿ = ( (1-p) e^-p¢u+pe^(1-p¢)u)ⁿ.

Now the calculus question: for what u ³ 0 the function f(u) = (1-p)e^-p¢u+pe^(1-p¢)u attains a minimum? The answer is u = 0 if p¢ £ p, or u = ln([(1-p)/ p][(p¢)/( 1-p¢)]). Therefore the minimal value is f(u_min) = [p/( p¢)]e^(1-p¢)u = exp(-p¢lnp¢/p-(1-p¢)ln(1-p¢)/(1-p)) ^[¯]
[ 2 If p = 1/2 then

Proof. This follows from the inequality (1+x)ln(1+x) +(1-x)ln(1-x) ³ x². ^[¯]

5.7  Conditional limit theorems

Suppose X_j are i.i.d. Conditional limit theorems say what is the conditional distribution of X₁, given the value of the empirical average ¹/_nå_{j = 1}ⁿ h(X_j). Such probabilities are difficult to simulate when ¹/_nå_{j = 1}ⁿ h(X_j) differs significantly from Eh(X).

5.8  Questions

Part 2
Stochastic processes

Chapter 6
Simulations

This chapter collects information about how to simulate special distributions.

Considerations of machine efficiency, in particular convenience of doing fixed precision arithmetics, make the uniform U(0,1) the fundamental building block of simulations. We do not consider in much detail how such sequences are produced - efficient methods are hardware-dependent. We also assume these are i. i. d., even though the typical random number generator returns the same pseudo-random sequence for each value of the seed, and often the sequence is periodic and correlated. This is again a question of speed and hardware only.

6.1  Generating random numbers

6.1.1  Random digits

The classical Peano curve actually maps the unit interval onto the unit square preserving the Lebesgue measure. Thus two independent U(0,1) are as ``random" as one!

Experiments with discrete outcomes aren't necessarily less random than continuous models. Expansion into binary fractions connects infinite tosses of a coin with a single uniform r. v.

Example 53 Let U be uniform U(0,2p). Random variables X_k = sign (sin(2^k U)) are i. i. d. symmetric independent.

Theorem 20 If x_j are independent identically distributed discrete random variables with values {0,1} and Pr(x = 1) = 1/2 then å_{k = 1}^¥ [1/( 2^k)]x_k in uniform U(0,1).

Proof. We show by induction that if U is independent of {x_j} uniform U(0,1) r. v., then [1/( 2ⁿ)]U+å_{k = 1}ⁿ [1/( 2^k)]x_k is uniform for all n ³ 0. For induction step, notice that in distribution

[1/( 2ⁿ)]U+å_{k = 1}ⁿ[1/( 2^k)]x_k @ ¹/₂x₁+¹/₂U. This reduces the proof to n = 1 case.

The rest of the proof is essentially the solution of Problem 2.18. Clearly Pr(¹/₂x₁+¹/₂U < x) = Pr(x₁ = 0)Pr(U/2 < x)+Pr(x₁ = 1)Pr(¹/₂+U/2 < x) Expression Pr(U/2 < x) is 0, 2x, or 1. Expression Pr(¹/₂+U/2 < x) is 0, 2x-1, or 1. Their average (check carefully ranges of x!) is

{

0	if x < 0
x	if 0 £ x £ 1
1	if x > 1

^[¯]
Coefficient 2 does not play any special role. The same fact holds true in any number system - if N > 1 is fixed, and U = åX_j N^-j is expanded in pase N, then X_j are independent uniform {0,1,..., N-1}-valued discrete random variables.

From the mathematical point of view all of the simulations can be based on a single U(0,1) random variable. In particular, to generate independently uniform integer numbers¹⁹ in the prescribed range 0 ... N we need to pick any u Î (0,1), and define X_j = N^ju mod N.

Clearly X_j is the integer part of NU_j, where U_j solve the recurrence

Many programs provide access to uniform numbers. These however might be platform-dependent, and are often of ``low quality". Often a person performing simulation may want to use the code they have explicit access to. What is ``random enough" for one application may not be random enough for another.

6.1.2  Overview of random number generators

There is good evidence, both theoretical (see (62)) and empirical, that the simple multiplicative congruential algorithm

The linear congruential method has the advantage of being very fast. It has the disadvantage that it is not free of sequential correlations between successive outputs. This shows up clearly in the fact that the consecutive k-points lie in at most N^1/k subspaces of dimension k-1.

Many system-supplied random number generators are linear congruential generators, which generate a sequence of integers n₁, n₂, ... each between 0 and N-1 by the recurrence relation

6.2  Simulating discrete r. v.

6.2.1  Generic Method - discrete case

This is rather easy to convert into the computer algorithm.

6.2.2  Geometric

6.2.3  Binomial

6.2.4  Poisson

An exact method to simulate Poisson distribution is based on the fact that it occurs the Poisson process, and that sojourn times in the Poisson process are exponential (see Theorem on page ). Therefore, to simulateX with Poiss(l) distribution, simulate independent exponential r. v. T₁, T₂,..., with parameter l = 1 and put as the value of X the first value of n such that T₁+...+T_n > l.

A reasonable approximation to Poisson Poiss(l) random variable X is obtained by simulating binomial Bin(n,p) random variable X¢ with l = np. Since X¢ £ n, use n large enough to exceed any realistic values of X. Run program LIMTHMS.BAS to compare the histograms - why is Poisson distribution more spread out than the binomial?

6.3  Simulating continuous r. v.

6.3.1  Generic Method - continuous case

Problem 42 Given a cumulative distribution function G(x) with inverse H(), and density g(x) = G¢(x), let U₁,U₂ be two independent uniform U(0,1) random variables. Show that X = H(U₁), Y = U₂g(X) are uniformly distributed in the region {(x,y): 0 < y < g(x)}.

6.3.2  Randomization

If the conditional density of Y given X = x is f(y|x), and the density of X is g(x), then the density of Y is òf(y|x)g(x) dx.

Example 54 Suppose Y has density cye^-y = Cå_{n = 0}^¥y(1-y)ⁿ/n!. Let Pr(X = n) = c/n!. Then the conditional distribution is f(y|X = n) = Cy(1-y)ⁿ, which is the distribution of a median from the sample of size 2n.

6.3.3  Simulating normal distribution

The exact simulation of normal distribution uses the fact that an independent pair X₁,X₂ of normal N(0,1) random variables can be written as

We simulate two independent normal N(0,1) r. v. from two independent uniform r. v. U₁, U₂ by taking Q = 2pU₁, R = Ö{-2lnU₂} and using formulas (65)-(66).

6.4  Rejection sampling

The name comes from the technique suggested by Problem 6.3.1. Instead of selecting points under the graph of f(x), we pick another function g(x) > f(x), which has known antiderivative with explicitly available inverse. We pick points under the graph of g(x), and "reject" those that didn't make it below the graph of f(x).

Often used density g(x) = c/(1+x²) leads to X = H(U₁), where H(u) = tan(pu/c).

Rejection sampling can be used to simulate continuous or discrete distributions. The idea behind using it in discrete case is to convert discrete distribution to a function of continous random variable. For example, to use rejection sampling for Poisson distribution, simulate the density f(x) = e^-ll^[x]/[x]! and take the integer part [X] of the resulting random variable.

6.5  Simulating discrete experiments

6.5.1  Random subsets

SUB RandomSubset( S())
n = UBOUND(S) ' read out the size of array
FOR j = 1 TO n
'select entries at random
S(j) = INT(RND(1) + 1)
NEXT j

For subsets of low dimension, a sequence of 0,1 can be identified with binary numbers. Set operations are then easily converted to binary operations on integers/long integers.

6.5.2  Random Permutations

Suppose we want to re-arrange elements a(1),..., a(n) into a random order. A quick method to accomplish this goal is to pick one element at a time, and set it aside. This can be easily implemented within the same sequence.

SUB Rearrange( A())
n = UBOUND(A) ' read out the size of array
ns = n 'initial size of randomization
FOR j = 1 TO n
'take a card at random and put it away
k = INT(RND(1) * ns + 1)
SWAP A(k), A(ns)
ns = ns - 1 'select from remaining a(j)
NEXT j

Problem 43 Let a₁,... a_n be numbers such that å_j a_j = 0, å_j a_j² = 1 Let X denote the sum of the first half of those numbers, after a random rearrangement.

• Find E(X), Var(X).
• Under suitable conditions, as n®¥, the distribution of X is asymptotically normal. Verify by simulations.
• Let Y be the sum of the first ¹/₄ of a_j after random rearrangement. Find E(X|Y) and E(Y|X).

6.6  Integration by simulations

Program listing on page pageref explains how to perform double integration by simulation. Here we concentrate on refining the procedure for single integrals by reducing the variance.

To evaluate J = ò_a^b f(x)  dx we may use the approximation ¹/_Nå_{j = 1}^N f(X_j)/g(X_j), where X_j are i.i.d. with the density g(x) such that ò_a^bg(x) dx = 1.

The error, as measured by the variance²¹, is

For smooth functions, a better approximation is obtained by selecting points more uniformly than the pure random choice. The so called Sobol sequences are based on sophisticated mathematics. Cookbook prescriptions can be found eg in Numerical recipes.

Note: The reliability of the Monte Carlo Method, and the associated error bounds depends on the quality of the random number generator. It is rather unwise to use an unknown random number generator in questions that require large number of randomizations. Example 55 The following problem is a challenge to any numerical integration method due to rapid oscillations, ò₀¹2sin²(1000 x)  dx = 1-[1/ 2000]sin2000 » 1.0004367

6.6.1  Stratified sampling

The idea of stratified sampling is to select different number of points from different subregions. As a simple example, suppose we want to integrate a smooth function over the interval [0,1] using n points. Instead of following with the standard Monte Carlo prescription, we can divide [0,1] into k non-overlapping segments and choose n_j points from the j-th subinterval I_j. An extreme case is to take k = n and n_j = 1 - this becomes a variant of the trapezoidal method. The optimal choice of n_j is to select them proportional to the local standard deviation of the usual Monte Carlo estimate of ò_Ij f(x)  dx. Indeed, denoting by F_j the estimator of the integral over I_j, the variance of the answer is Var(å_{j = 1}^k F_j) = å_j Var(F_j) = å_j s_j²/n_j. The minimum under the constraint å_jn_j = n is n_j » s_j.

A simple variant of recursive stratified sampling is to generate points and subdivisions based on estimated values of the variances.

6.7  Monte Carlo estimation of small probabilities

Example

Pr(X₁+...+X_n > an) = ò_{y1+...+yn > an}e^-ly₁... e^-ly_nf(y₁)... f(y_n) dy₁... dy_n. This leads to the following procedure.

• Simulate N independent realizations of the sequence Y₁,..., Y_n.
• Discard those that do not satisfy the constraint y₁+...+y_n > an.
• Average the expression e^-lY₁... e^-lY_n over the remaining realizations to get the desired estimate.

Example 56 Suppose X_j are {0,1}-valued so that X₁+... X_n is binomial Bin(n,p). What is the distribution of Y_j? Simplify the expression e^-lY₁... e^-lY_n.

[ 26 Write the program computing by simulation the probability that in a n = 10 tosses of a fair coin, at least 8 heads occur. Once you have a program that does for n = 10 a comparable job to the ``naive" program below, try Exercise 1.5 on page pageref.

Here is a naive program, that does the job for n = 10, but not for n = 100. It should be used to test the more complex ``tilted density" simulator.

PROGRAM heads.bas
'

'Simulating N fair coins
' declarations
DECLARE FUNCTION NumHeads% (p!, n%)
DEFINT I-N  ' declare integer variables

'prepare screen
CLS

PRINT "Simulating toss of n fair coins"

'get users input
n = 10 'number of trials 
INPUT "Number of coins n=", n
pr = .5 ' fairness of a coin
frac = .8 ' percentage of heads seeked
' get the frac from the user
PRINT " How often we get more than f heads? (where 0<f<"; n; ")"
INPUT "f=", frac
IF frac >= 1 THEN frac = frac / n 'rescale if too large
'tell user what is going on
PRINT "Hit any key to see the final answer"
LOCATE 20, 10
PRINT "With some patience you may see digits stabilize"

DO 'main loop
        T = T + 1
        IF NumHeads(pr, n) > frac * n THEN s = s + 1
        IF INKEY$ > "" THEN EXIT DO
        LOCATE 10, 10
        PRINT "Trial #"; T
        PRINT "Current estimate"; USING "##.#####"; s / T
LOOP
'print the answer
PRINT
PRINT "Prob of more then "; INT(frac * n); " heads in "; n; " trials is about "; s / T
END

DEFINT A-H, O-Z
FUNCTION NumHeads (p!, n)
'simulate a run of n coins
 h = 0
FOR k = 1 TO n
        IF RND(1) < p! THEN h = h + 1
NEXT k
NumHeads = h
'

Chapter 7
Introduction to stochastic processes

stochastic, a. conjectural; able to conjecture
Webster's New Universal Unabridged Dictionary

Stochastic processes model evolution of systems that either exhibit inherent randomness, or operate in an unpredictable environment. This unpredictability may have more than one form, see Section .

Probability provides models for analyzing random or unpredictable outcomes. The main new ingredient in stochastic processes is the explicit role of time. A stochastic process is described by its position X(t) at time t Î [0,1], t Î [0,¥), or t Î {0,1,...}.

From the conceptual point of view, stochastic processes that use discrete moments of time t Î {0,1,...} are the simplest. Since the discrete moments of time can represent arbitrarily small time increments, discrete time models are rich enough to model real-world phenomena. Continuous time versions are convenient mathematical idealizations.

7.1  Difference Equations

Mathematical models of time evolution of deterministic systems often involve differential equations.

Difference equations are discrete analogues of the differential equations. Difference equations occur in applied problems, and also when solving differential equations by series expansions, or by Euler's method. The concepts of numerical versus analytical solution, initial values, boundary values, linearity, and superposition of solutions occur in the discrete setup in complete analogy with the theory of differential equations.

A difference equation determines an unknown sequence (y_n) through a relation that specifies the pattern. The relation often has the form of an equation that ties the consecutive values of the sequence. The technical term for such an equation is recurrence.

We will consider only special cases of classes of equation

7.1.1  Examples

Example 57

Suppose a sequence y_n is to satisfy

It is easy to write a short program that computes the values of y_n. But to compute the actual sequence we need to specify the initial value y₀. Table gives sample outputs of such a program for several choices of y₀.

y0	y1	y2	y3	y4	y5	y6	y7	y8
-1	.5403023	.8575532	.6542898	.7934803	.7013688	.7639596	.7221025	.7504177
-.5	.8775826	.6390125	.8026851	.694778	.7681958	.7191654	.7523558	.7300811
0	1	.5403023	.8575532	.6542898	.7934803	.7013688	.7639596	.7221025
.5	.8775826	.6390125	.8026851	.694778	.7681958	.7191654	.7523558	.7300811
1	.5403023	.8575532	.6542898	.7934803	.7013688	.7639596	.7221025	.7504177
1.5	.0707372	.9974992	.542405	.8564697	.6551088	.7929816	.7017242	.7637303
2	-.4161468	.9146533	.6100653	.8196106	.6825058	.7759946	.7137247	.7559287

Similar numerical procedures show up in differential equations, where they are used to approximate continuous solutions by discretizing time. The procedure is called Euler method, or tangent line method.

Example

On the other hand, we may notice that equation (68) defines the arithmetic progression. Instead of a table like Table 7.1, we can write down the solution for all possible values of y₀ and for all n. Namely,

Example 58 Here is another well known difference equation with obvious solution. Suppose

Examples 7.1.1 and 7.1.1 are deceptively simple. Not every difference equation has a simple or easy to guess answer.

Example

7.2  Linear difference equations

Many interesting difference equations, including (68), (70), (72) fall into the category of linear difference equations with constant coefficients. The first order linear difference equations with constant coefficients has the form

A method to solve difference equation of order 2 consists of the following steps:

• First we solve the ``homogeneous equation" y_n+2+ay_n+1+b y_n = 0. This is accomplished through the substitution y_n = rⁿ,which leads to the characteristic equation r²+ar+b = 0 for r. Once two roots are found, the general solution is y_n = C₁r₁ⁿ+C₂r₂ⁿ. The notable exception when the roots are equal is y_n = C₁r₀ⁿ+C₂nr₀ⁿ.

• Secondly, we find any solution of the non-homogeneous equation, disregarding initial conditions. This can be accomplished by the method of varying parameters: try y_n = C₁(n)r₁ⁿ+C₂(n)r₂ⁿ. Or by guessing. In the often encountered polynomial case g(x) = n^l the trial solution y_n = n^s(u₀+u₁n+...+u_l n^l) will work (s = 0 except when y_n = n solves the homogeneous equation, in which case s = 1).
• In the final step the general solution to homogeneous equation is combined with the particular solution to non-homogeneous equation, and the initial value problem is solved.

Example

Suppose you borrow y₀ dollars on fixed monthly interest rate r. If you do not make any payments on your loan, then your balance will ``balloon" exponentially. Formula y_n = (1+r)ⁿy₀ expresses monthly balance after n months when no loan payments are made. Some people prefer to pay a fixed amount of p dollars at the end of each month. If p is large enough, then the balance may even shrink down! This situations is easily described by a difference equation. To write it down, compute the next month balance, if the previous month balance is known:

In differential equations, we often solve a similar problem in continuous time, with continuous compounding and continuous payments schedule. This is a mathematical simplification of the actual banking situation, but the answers are reasonably close, and they are easier to get.

Here is an example of the mathematical problem that leads to a natural, but non-trivial²⁶ difference equation.

Example 59 Suppose we want to find the formula for the sum of all consecutive integers from 1 to n. Let the answer by y_n. Then the recurrence formula

Now (74) can be written as y_n+1-y_n = n+1 and actually this is why we use the name difference equations rather than recurrence relations. Notice that we do know the solution of its continuous analogue. Equation y¢ = t+1 for unknown function y = y(t) resembles (74). Its solution²⁷ is y(t) = ¹/₂t²+t. So to find the answer to (74) we may try substituting u_n = ¹/₂n²+n into the equation. Unfortunately, simple arithmetics shows that we don't get the right answer, as

The method we used - subtracting the equations - works for the so called linear difference equations (and also linear differential) equations. It is closely related to the Principle of Superposition for linear differential equations.

7.2.1  Problems

1. Check if the given sequence solves the difference equation.

   1. Equation: y_n+1 = -y_n. Sequence y_n = cosnp.
   2. Equation: y_n+1 = y_n+2y_n-1. Sequence y_n = 2ⁿ.
   3. Equation: y_n+1 = y_n+2y_n-1. Sequence y_n = -¹/₂n³+3n²-⁵/₂n+1.
   4. Equation: y_n+1 = 2y_n-y_n-1+2. Sequence y_n = n².

2. Write down the difference equation that you need to solve each of the following problems (do not solve the equation, nor the problem).

   1. A loan of $1000 has interest rate that varies with time as follows: first month interest is 0%, second month interest is [1/ 12]%, third month interest is [2/ 12]%, etc with monthly interest increasing by [1/ 12]% every month. Determine the fixed monthly payment p that will pay this loan within one year.

   2. A loan of $1000 has constant monthly interest rate of [1/ 12]%. I arranged my monthly payments in the following fashion: at the end of the first month I will pay nothing, at the end of the second month I will pay $10, at the end of the third month I will pay $20, etc, with monthly payments increasing by $10 every month. Determine when I will pay off this loan and how much money I will pay in total.

3. Find the general solution of the following difference equations.

   1. y_n+1 = ¹/_n y_n
   2. y_n+1 = [n/( n+1)]y_n
   3. n(n+1)y_n+1 = y_n

4. Solve the given initial value problem for the difference equation.

   1. y_n+1 = ¹/_n y_n; y₀=1
   2. y_n+1 = [n/( n+1)]y_n-1 ; y₀ = 1, y₁ = 0
   3. n(n+1)y_n+1 = 2y_n ;y₀ = 1, y₁ = 0

5. Find the formula for

   1. y_n = 1²+2²+...+n²
   2. y_n = 1³+2³+...+n³
   3. y_n = ¹/₂+¹/₆+...+[1/( n(n+1))]
   4. y_n = 1+r+r²+...+r^n-1
   5. y_n = r+2r²+3r³+...+nrⁿ

6. Each solution of equation (67) has the limit lim_n®¥y_n. Show that the first digits of this limit are .739085133.

7.3  Recursive equations, chaos, randomness

Before jumping to models that use explicit randomness in their evolution, it is quite illuminating to analyze first some mathematically simple and well defined ``deterministic evolutionary processes". The following set of examples describes deterministic evolution in discrete time of a system described by a single numerical parameter x Î (0,1). All of the examples fall into the category of (non-linear) difference equations: given initial value x₀ Î (0,1) and a simple evolution equation of the form x_n+1 = g(x_n), we are supposed to make inferences about the behavior of the solutions.

There is no randomness in the evolution itself. But since we are allowed to choose any initial condition, and no initial condition can be measured exactly, we may as well consider the initial value to be random. Example 60 Equation

Example 61 Let {x} denote the fractional part of x. Equation

Here is a rather surprising fact. Suppose x₀ is selected at random. Since x_k is a function of x₀, and x₀ is random, x_k becomes a random variable. It may happen that x_k < 1/2. Call this event A_k. The reasoning presented in Section 6.1.1 implies that events {A_k} are independent and have the same probability Pr(A_k) = 1/2. Therefore the deterministic evolution equation contains at least as much randomness as the tosses of a coin.

The last example may perhaps give the impression that it is the discontinuity of the evolution equation that is the source of difficulties. This is not at all the case.

Example 62 Equation

Example 7.3 indicates that ``naive" prediction of the future of a system is unreliable even within the realm of deterministic evolution equations.

On the other hand, there are aspects of the evolution that we can analyze reliably. These deal with the average behavior of the evolution.

There are two classes of questions that we can answer, but both deal with statistical nature of the evolution:

• what happens if the same experiment is repeated many times (with slightly different initial conditions)?
• What is the average behavior of the system over long periods of time?

• What is the average ¹/_nå_{j = 1}ⁿx_j?
• What is ¹/_nå_{j = 1}ⁿU(x_j) for a given function U?
• How often x_j < 1/2? (Meaning - what proportion of j £ T satisfies the condition for large T.)

• How often x₃ < 1/2 when x₀ is selected according to density g(x)? (Meaning - what proportion of x₀ satisfies the condition for large number of initial points x₀.)

7.4  Modeling and simulation

The quick way to get insights into operation of real systems is to model their behavior. Here are examples of enterprises that operate under randomness. Mathematicians devised methods of modelling each of these. But it is interesting also to simulate their behavior. Notice that simulations require assumptions, but so do the analytical methods. Regardless of the method, we have to be careful about what are the questions we can answer.

Example 63 Suppose we want to study how much capital is needed to run a casino. We need to answer the following.

• How often do people win?
• How likely is the casino to loose money in a day? In a month?
• How much capital should be kept on hand to cover the losses?

Example 64 Suppose we want to study how much capital is needed to insure cars. We need to answer the following.

• How often do accidents occur? How expensive are repairs/medical costs?
• How likely is that the insurance company looses money in a day? In a month?
• How much capital should be kept on hand to cover losses?

There is also a number of questions that we do not want to answer. For instance we do not want to predict whether I'll file an insurance claim today, driving back home on I-74 without enough sleep since I was preparing this class until 3AM.

Example 65 A construction company has n jobs to be performed in the future. For each of these jobs, experts provide estimate of the cost. Then, after questioning, they complement the estimates by the lower/upper bound for the costs. How much money should be allotted?

Example 66

A store averages l(t) customers per hour at time t of the day. Each customer brings some (random) profit. However, a customer may just leave the store without shopping, if the lines are too long.

• How many cashiers should be available for each shift (time) t?
• What are the profits on average?

7.5  Random walks

A random walk is a process of the form X_n = å_{j = 1}ⁿ x_j, where x_j are i. i. d. Random walks have independent increments , and describe the accumulation of independent contributions over time.

Random walks are examples of Markov processes which will be studied in detail in Chapter . Their special structure allows to analyze them independently of the general theory.

7.5.1  Stopping times

Definition 19 T:W® INÈ¥ is a stopping time, if the event {T = n} is independent of x_n+1,x_n+2,....

This definition is specialized to random walks. In a more general Markov case the definition is less transparent but captures the same idea.

The most important example of a stopping time is the first entrance time T = inf{k: X_k Î A}. An example that is not a stopping time is the last exit from a set.

When T < ¥ we define random sums X_T = å_{j £ T}x_j. The following theorem is an exercise when T is independent of x.

Theorem 21 If x_j are i. i. d., Ex = m, and ET < ¥ is a stopping time then

Proof. EX_T = å_n EX_n Pr(T = n) = å_n å_{k = 1}ⁿEx_k Pr(T ³ k). Since x_k and {T ³ k} = {T < k}¢ are independent, therefore ES_T = må_n Pr(T ³ n) = mET by tail integration formula (21). ^[¯]
Theorem 22 If x_j are i.i.d., Ex = m, Var(x) = s² < ¥, and ET < ¥ then

(These formulas are of interest in branching processes, and in chromatography.)

Example 67 The number of checks cashed at a bank per day is Poisson random variable N with mean l = 200. The amount of each check is a random variable with a mean of $30 and a standard deviation of $5. If the bank has $6860 on hand, is the demand likely to be met?

Problem 44 Suppose x_k, T are independent. Find the variance of X_T in terms of the first two moments of x, T.

7.5.2  Example: chromatography

The sample is injected into a column, and the molecules are transported along the length by electric potential, flow of gas, or liquid. The basis for chromatographic separation of a sample of molecules is difference in their physical characteristics. The molecules switch between two phases: mobile, and stationary, and the separation of compounds is caused by the difference in times spend in each of the phases.

Suppose that the molecules of a compound spend random independent amounts of time U_k in mobile phase and random amount of time W_k in the stationary phase. Thus at time t the position of a molecule is given by a random sum vå_{j = 1}^T(t) U_j, where T(t) = inf{k: å_{j = 1}^kU_j+W_j > t}.

Section 7.5.1 gives formulas for the mean and the variance of the position. Since the number of transitions T is likely to be large for a typical molecule, it isn't surprising that the actual position has (asymptotically) normal distribution. (The actual Central Limit Theorem for random sums is not stated in these notes.)

[ 27 Simulate the output of the chromatography column of fixed length separating a pair of substances that have different distributions of mobile and stationary phases U_k, W_k. Select a hundred particles of each substance, and measure the degree of separation at the end of the column.

7.5.3  Ruining a gambler

Suppose a gambler can afford to loose amount L > 0, while the casino has capital C < 0. Let x_j = ±1 be i. i. d. random variables modelling the outcomes of consecutive games, S_n be the partial sums (representing gains of the gambler), and let T = inf{k:S_k ³ L or S_k £ C} be the total number of games played. Then Pr(T > k) = Pr(C < S_k < L) and thus ET = å_kPr(C < S_k < L).

The special case Pr(x = ±1) = 1/2 is easily solved, since in this case ES_T = ExET = 0. Let p = Pr(S_T = C) denote the probability of ruining the casino. Since S_T is either C, or L we have 0 = ES_T = pC+(1-p)L, giving p = L/(L-C). This formula means that a gambler has a fair chance of ruining a casino in a fair game, provided he brings with him enough cash L.

For more general random walks (and less fair games) probability of gambler's ruin can be found explicitly using the one-step-analysis (Section ). It is also interesting to find how long a game like that would last on average. (The expression for ET given above is not explicit.)

7.5.4  Random growth model

Let X_t denote the number of bacteria in t-th generation, with X₀ = 1. Assume that a bacteria can die with probability q > 0, or divide into two cells with probability p = 1-q, and that all deaths occur independently. Our goal here is to find the average number of bacteria m(t) = E(X_t) in the t-th generation. This can be recovered from Theorem 7.5.1. Instead, we show another method based on conditioning.

The number of bacteria in the next generation is determined by binomial probabilities: Pr(X_t+1 = 2k|X_t = n) = (ⁿ_k)p^k q^n-k. Therefore E(X_t+1|X_t) = 2p X_t and the average population size m(t) = E(X_t) satisfies difference equation m(t+1) = 2p m(t). We have m(t) = (2p)^t. In particular, the population on average grows exponentially when p > 1/2.

It is perhaps surprising that a population of bacteria with p = 3/4, which on average grows by 50% per generation, has still a ¹/₃ chance of going extinct. One way to interpret this is to say that infections by a ``deadly" and rapidly developing desease may still have a large survival rate without any intervention of medicine, or immune system. (The methods to compute such probabilities will be introduced in Section . The answer above assumes infection by a single cell.)

Problem 45 Find the formula for the variance Var(X_t) of the number of bacteria in t-th generation.

Problem 46 What is the probability that an infection by 10 identical bacteria with the doubling probability p = 3/4 dies out?

Chapter 8
Markov processes

evolution, n. [L. evolutio (-onis)], an unrolling or opening...
Webster's New Universal Unabridged Dictionary

Markov processes are perhaps the simplest model of a random evolution without long-term memory.

Markov process is a sequence X_t of random variables indexed by discrete time t Î ZZ_{\scriptscriptstyle +} , or continuous t ³ 0 that satisfies the so called Markov property. The set of all possible values of variables X_t is called the state space of the Markov chain. Typical examples of state spaces are IR, IN, the set of all non-negative pairs of integers, and finite sets.

Markov chains are Markov processes with discrete time. Thus a Markov chain is an infinite sequence {X_k}_{k Î ZZ\scriptscriptstyle +} of (usually, dependent) random variables with short-term (one-step) memory.

8.1  Markov chains

The formal definition of the Markov property is as follows. Definition 20 A family of discrete r. v. {X_k}_{k Î ZZ\scriptscriptstyle +} is a Markov chain, if

Examples of Markov chains are:

• A sequence of independent r. v.
• A constant random sequence X_k = x.
• Random walks (sums of independent random variables).

Examples of non-Markov processes are easy to construct, but lack of Markov propertry is not obvious to verify. In general, if X_k is a Markov process, Y_k = f(X_k) may fail to be Markov.

8.1.1  Finite state space

Markov condition

Example 68 Suppose X_n is a Markov chain with periodic transition probabilities P_n = P_n+T. Then Y_n = (X_n+1,X_n+2,..., X_n+T) is a homogeneous Markov chain.

Problem 47 Suppose x_j are independent {0,1}-valued with Pr(x = 1) = p. Let X_n = ax_n+bx_n+1, where ab ¹ 0.

Explain why X_n is a Markov chain.

Write the transition matrix for the Markov chain X_n.

[ 3 The probabilities Pr(X_n = j|X₀ = i) are given by the i,j-entries of the matrix Pⁿ

Proof. This is the consequence of Markov property (78) and the total probability formula (10). ^[¯]
Powers of moderately sized matrices are easy to compute on the computer. Section indicates a mathematical method of computing Pⁿ for small dimensions using the Cayley-Hamilton theorem. Under certain conditions the powers converge. [ 28 Find lim_n®¥Pⁿ for the matrix from Problem 8.1.1.

8.1.1 Stationary Markov processes