AMTH142 Lecture 13. Random Numbers - PDF

Please download to get full document.

View again

of 15
All materials on our website are shared by users. If you have any questions about copyright issues, please report us to resolve them. We are always happy to assist you.
Information Report

Health & Lifestyle


Views: 0 | Pages: 15

Extension: PDF | Download: 0

Related documents
AMTH142 Lecture 13 Random Numbers April 20, 2007 All the Scilab functions defined in this Lecture can be found in the file l13.sci in the directory for this lecture. Contents 13.1 Review of Probability
AMTH142 Lecture 13 Random Numbers April 20, 2007 All the Scilab functions defined in this Lecture can be found in the file l13.sci in the directory for this lecture. Contents 13.1 Review of Probability Definitions Uniform Density Normal Density Random Number Generators Algorithms Scilab Some Common Distributions Uniform Distribution Normal Distribution Uniform Discrete Distributions Non-Uniform Random Numbers Inverse Transform Method Acceptance/Rejection Method Special Distributions 13.1 Review of Probability This is just a reminder of few facts about continuous probability distributions covered in MATH Definitions Continuous probability distributions are described by probability densities, i.e. real valued functions ρ(x) with the properties: 1. ρ(x) ρ(x)dx = 1. by 3. Prob(a x b) = b a ρ(x)dx. The distribution function associated with the density ρ(x) is defined It has the properties: 1. F (s) = Prob(x s). 2. lim x F (x) = 0 3. lim x F (x) = 1 F (x) = x 4. F (x) is an increasing function of x. 5. df dx = ρ(x) It follows from 1. that and ρ(t)dt. Prob(a x b) = F (b) F (a). The mean and variance of a density ρ(x) are defined by µ = Mean(x) = σ 2 = Var(x) = xρ(x)dx (x µ) 2 ρ(x)dx The standard deviation, σ, is the square root of the variance. 2 Uniform Density The uniform density on the interval [a, b] is defined by 0 x a ρ(x) = 1 b a a x b 0 x b An important special case is the uniform density on [0, 1] 0 x 0 ρ(x) = 1 0 x 1 0 x 1 This density has mean µ = 1/2 and variance σ 2 = 1/ Normal Density The normal density with mean µ and standard deviation σ is defined by ρ(x) = 1 σ (x µ) 2 2π e 2σ 2 The case µ = 0, σ = 1 is especially important ρ(x) = 1 e x2 2 2π 13.2 Random Number Generators Computer algorithms for generating random numbers are deterministic algorithms. Although the sequence of numbers produced by a random number generator appears random, the sequence of numbers is completely predictable and for this reason they are often called pseudo-random. Since computers have only a finite number of different states, the sequence of pseudo-random numbers produced by any deterministic algorithm must necessarily repeat itself. The number of random numbers generated before the sequence repeats is called the period of the generator. A good random number generator should have the following characteristics: 1. Randomness. It should pass statistical tests of randomness. 2. Long Period. For obvious reasons. 3. Efficiency. This is important since simulations often require millions of random numbers. 3 4. Repeatability. It should produce the same sequence of numbers if started in the same state. This allows the repetition and checking of simulations. Almost all random number generators used in practice produce uniform [0, 1] distributed random numbers and from these random numbers with other distributions can be produced if required Algorithms There are many algorithms for generating pseudo-random numbers. Linear Congruential Generators The simplest are the linear congruential generators. Starting with the seed x 0, these generate a sequence of integers by x k = (ax k 1 + c) mod M (1) where a, c and M are given integers. All the x k are integers between 0 and M 1. In order to produce floating point numbers these are divided by M to give a floating point number in the interval [0, 1). Below is a Scilab function 1 implementing a linear congruential generator. Here n is the number of terms generated a, c and m are as in equation (1) and x0 is the initial value or seed. A vector of n+1 values including the seed is returned. function x = lcg(n, a, c, m, x0) x = zeros(1,n+1) x(1) = x0 for i = 2:n+1 x(i) = pmodulo(a*x(i-1)+c, m) end x = x/m endfunction The quality of a linear congruential generator depends on the choice of a, b and M, but in any case the period of such a generator is at most M (why?). As we will see later, Scilab s basic rand generator is a linear congruential generator. The following example illustrates a common problem with some random number generators. We will take the linear congruential generator with a = 1203, c = 0, m = With initial value x 0 = 1, this generator has period 512. We will run through a full period of the generator: 1 The version in l5.sci has been modified to avoid a problem with rounding error. 4 -- xx = lcg(511, 1203, 0, 2048, 1); -- xx(1:10) ans = column 1 to 4! ! column 5 to 8! ! column 9 to 10! ! Now take successive pairs of values as the x and y coordinates of a point in the plane and plot the results. -- x = xx(1:2:511); -- y = xx(2:2:512); -- plot2d(x,y,style = -1) This latticing is a common problem with random number generators. Of course, we have used a very simple generator for which the problem is obvious, for other generators the problem may occur in higher dimensions. That is, taking successive n-tuples of values and plotting them (hypothetically at least) in n-dimensional space the points may tend to lie on a lattice of lower dimensional subspaces. 5 Fibonacci Generators Another common type of generator are the Fibonacci generators. typical example generates a sequence by A x k = (x k 17 x k 5 ) mod 1 This example uses the previous 17 numbers (and needs 17 seeds to get started). It has a number of advantages compared to congruential generators: 1. The x k are floating point numbers, so no division is needed. 2. It has a much longer period since the repetition of a single number does not entail repetition of the whole sequence. 3. It can have much better statistical properties. Combined Random Number Generators By combining a few simple random number generators we can obtain a generator with better properties. One generator used by Scilab combines four linear congruential generators t n =45991t n 1 mod u n =207707u n 1 mod v n =138556v n 1 mod w n =49689w n 1 mod x n =(t n u n + v n w n ) mod This generator has a period of about Other Random Number Generators Modern random number generators are based on number theory and have quite sophisticated implementations. One in common use is the Mersenne twister which takes a vector of 625 as its seed Scilab Scilab has two random number generators rand and grand. We will only discuss the simpler generator rand which can produce either uniform [0, 1] or normal µ = 0, σ = 1 random numbers. Uniform random numbers are the default. We have already seen how rand works: 1. rand(m, n) gives a m n matrix of random numbers. 6 2. rand(m, n, normal ) gives a m n matrix of normally distributed random numbers. 3. rand(a) or rand(a, normal ) for matrix a gives a matrix of random numbers the same size as a. 4. rand( seed, 0) resets the random number generator to its original state. This is handy if you want to repeat an experiment using the same random numbers. As mentioned earlier, rand is a linear congruential generator: x n = (ax n 1 + c) mod M with a = , c = , and M = Let us check this. First we set the seed of rand to 0 (which is the seed on the first call to rand) and generate terms. -- rand( seed ,0) -- x1 = rand(1,10000); Now generate terms of lcg 2 with parameters above and seed 0. -- x2 = lcg(10001, , , 2^31, 0); We drop off the initial term (which is the seed) and look at first few terms: -- x2 = x2(2:10001); -- x1(1:8) ans = column 1 to 4! ! column 5 to 8! ! -- x2(1:8) ans = 2 Use the version in l5.sci, see previous footnote. 7 column 1 to 4! ! column 5 to 8! ! and finally find how many terms differ: -- sum(x1 ~= x2) ans = Some Common Distributions Two simple change of variable tricks are very important Uniform Distribution If x is a uniform [0, 1] distribution, then the change of variable y = (b a)x + a gives a uniform [a, b] distribution. Here is an example generating a uniform [ 5, 5] distribution: -- x =rand(1,10000); -- y = 10*x - 5; -- histplot(100,y) Normal Distribution If x has a normal µ = 0, σ = 1 distribution, then the change of variable y = ax + b gives a normal distribution with µ = b and σ = a. Here is an example generating a normal µ = 100, σ = 10 distribution: -- x = rand(1,10000, normal ); -- y = 10*x + 100; -- histplot(100,y) Uniform Discrete Distributions It is quite common to want to generate uniformly distributed random integers. Typically they will be in the range 0 k n or 1 k n. These can be obtained from uniform [0, 1] random floating numbers x by 1. floor((n+1)*x) gives integers in the range 0 k n. 2. floor(n*x+1) gives integers in the range 1 k n. Here is an example generating integers in the range 1 to 10 (inclusive): -- x = rand(1,30); -- y = floor(10*x+1) y = column 1 to 10 9 ! ! column 11 to 20! ! column 21 to 30! ! 13.4 Non-Uniform Random Numbers Inverse Transform Method Consider a probability density ρ(x) and let F (x) be the corresponding distribution function. Then F : R [0, 1] and since F (x) is an increasing function it has an inverse function G(y) with G : [0, 1] R x G(y) b a F (a) F (b) 1 y Let y be uniformly distributed on [0, 1] and set x = G(y). Then Prob(a x b) = Prob(F (a) y F (b)) = F (b) F (a) 10 Where the last equality follows from the uniformity of y. The statement Prob(a x b) = F (b) F (a) is equivalent to the statement that x has distribution function F (x) and hence x has density ρ(x). In summary: Given a probability density ρ(x) with distribution function F (x), let G(y) be the inverse function of F (x). Then if we take y to be uniformly distributed on [0, 1], x = G(y) has probability density ρ(x). Exponential Distribution This is defined by ρ(x) = { 0 x 0 λe λx x 0 with distribution function F (x) = { 0 x 0 1 e λx x 0 (Here λ is a parameter). Setting y = F (x) = 1 e λx and solving for x we obtain the inverse function x = G(y) = ln(1 y). λ Lets see how it works in Scilab; we will take λ = 2. -- y = rand(1,10000); -- x = - log(1-y)/2; -- histplot(100,x) We can compare the histogram to exact density: -- xx = 0:0.01:5; -- plot2d(xx, 2*exp(-2*xx)) 11 One drawback of the inverse transform method is that we need to know the inverse of the distribution function. For most distributions there is no explicit formula for this inverse and we must resort to approximations or seek another method Acceptance/Rejection Method This method too has its drawbacks. We will need to assume that we are working with a density which is zero outside of a finite interval. Suppose that ρ(x) is zero outside of the interval [a, b] and furthermore that ρ(x) is bounded above by c. Reject y c Accept ρ(x) a b x Now generate points (x i, y i ) with x i uniformly distributed in [a, b] and y i uniformly distributed in [0, c]. If y i ρ(x i ) then we accept the value x i, if y i ρ(x i ) then we reject the value x i. The values x i which are accepted will have the probability density ρ(x). 12 An Example Consider the hat-shaped probability density defined on [ 1, 1] by { x x 0 ρ(x) = 1 x 0 x 1 1 ρ(x) 1 1 We want to write a Scilab function to generate n random numbers with this density using the acceptance/rejection method. One thing to keep in mind when using the acceptance/rejection is that we do not know before hand how many random numbers we need to generate. First we need the density function 3 : function y = rhohat(x) if (x 0) then y = x + 1 else y = 1 - x end endfunction Now the random number generator: function x = randhat(n) x = zeros(1,n) k = 0 // keep count of numbers generated while (k n) xx = *rand(1,1) // uniform on [-1,1] yy = rand(1,1) if (yy = rhohat(xx)) // accept xx k = k + 1 x(k) = xx end end endfunction 3 The version of rhohat in l5.sci has been modified to allow a vector x as the argument. 13 And testing it -- x = randhat(10000); -- histplot(100, x) -- xx = -1:0.01:1; -- plot2d(xx, rhohat(xx)) Special Distributions The two methods discussed above, acceptance/rejection and inverse transform, when applied with some ingenuity can be used to generate random numbers with almost any distribution. However for most distributions there are specific techniques that perform better than these more general techniques. The Scilab function grand has facilities for generating random numbers from most of common distributions. Normal Distribution Here is a clever method for generate normally distributed random numbers: generate two uniform [0, 1] random numbers x 1 and x 2, then y 1 = sin(2πx 1 ) 2 ln x 2 and y 2 = cos(2πx 1 ) 2 ln x 2 are both normally distributed with µ = 0, σ = 1. 14 Here is how it works in Scilab: function x = randnorm(n) // assumes n even m = n/2 x1 = rand(1,m) x2 = rand(1,m) y1 = sin(2*%pi*x1).*sqrt(-2*log(x2)); y2 = cos(2*%pi*x1).*sqrt(-2*log(x2)); x = [y1 y2] endfunction -- x = randnorm(10000); -- histplot(100,x) We can compare the histogram to exact density: -- xx = -4:0.01:4; -- yy = exp(-(xx.^2)/2)/sqrt(2*%pi); -- plot2d(xx, yy)
We Need Your Support
Thank you for visiting our website and your interest in our free products and services. We are nonprofit website to share and download documents. To the running of this website, we need your help to support us.

Thanks to everyone for your continued support.

No, Thanks