Project one. Problem Three. Mathematics 502 Probability and Statistics
Nasser Abbasi, September 26,2007. California State University, Fullerton
Problem 3 part (a)
We are asked to generate R.V's from . We note as shown in the problem itself, that R.V. X can be written as product of 2 R.V WT where Wis ±1 with probability each. Hence to generate R.V. we do the following. We generate n R.V. from uniform distribution [0,1] using Mathematica random number generator. Then we check if each number is or not, and we generate 1 or -1 as the case may be. We then generate n random variables from the exponential distribution, which we know how to do from part (a). Then we multiply the above 2 vectors, element wise, with each others. The first vector being the vector of 1's and -1's. And the second vector being the RV's from the exponential distribution. This is the algorithm
Algorithm
Input: λ,n (number of random variables to generate)
output: list of random numbers which belong to density
Seed the random number generator with unique value for us.
A = Generate n random numbers from the exponential distribution with parameter λ (CALL problem 1 part(a) with the input λ,n) This uses method and uniform random number generator as well.
B = Generate n random numbers from uniform random number generator [0,1]
FOR i in 1..n LOOP -- Note: This is algorithm view. In code 'vectored' operation is used.
IF B(i)<.5 THEN
B(i) = 1
ELSE
B(i) = -1
END IF
END LOOP
result = B * A
Now generate a histogram from the result above.
The following function implements the above algorithm
Code Implementation
Define the function which was derived earlier. This is the inverse of the CDF of the exponential density function
In[208]:=
In[210]:=
In[211]:=
In[212]:=
Test the above function by plotting the histogram generated for say n=10000 overlaid by the true double exponential density function.
First, define the double exponential function
In[213]:=
Now do the overlay plot
This function makes a histogram which is scaled to be used to overlay density plots, or other functions.
Input: originalData: this is an array of numbers which represents the data to bin
nBins: number of bins
output: the histogram itself but scaled such that area is ONE
In[214]:=
In[216]:=
Out[219]=
Problem 3 part(b)
In this part, we need to generate a list of r.v's that belong to normal distribution N(0,1), using uniform random number generator U[0,1] and using the random numbers generated from the double exponential density function in part (a) above. We are asked to use the accept/reject method.
First the method is explained, then the algorithm outlined, then the implementation shown and a test case given, then a GUI interface written to test the algorithm for different parameters values.
Accept/Reject algorithm
input: n (number of random variables to generate)
λ (the exponential density parameter)
f(x) the density function for random variable X which we wish to generate random variables
(x) the density which we will use to help in generating the random variables from (x). This density is such that it is easy to generate random variables from. Much easier that from f(x) and that is why it was selected.
output: list of random numbers of length n from f(x)
Step 1: Find C. Where To solve this, this is the algorithm
Algorithm for step 1: Let (since double exponential is symmetric, I'll use one sided version). Let . Now find the ratio Now find where the maximum of this ratio is using normal calculus method: Take the derivative w.r.t. x and set it to zero. Solve the resulting equation for x. Evaluate the ratio at this root. This gives C. We find that C=1.31549 The following few lines of code finds C:
In[220]:=
Out[223]=
Step 2: Now that we found C in step 1, then the envelop function becomes
Step 3: seed the number random generator
initialize an array d of size n to contain all the accepted random numbers generated
initialize counter number_accepted=0
WHILE number_accepted < n DO
generate r.v. from U[0,1] call it u.
Generate r.v. from double exponential density (using part(a)) call this x
IF u*C* THEN
d[i]=x
number_accepted++
END IF
END LOOP
Step 4: Now array d contains the n random numbers generated from the normal density N[0,1]. Make histogram and overlay it over N[0,1]
Diagram showing main steps in the algorithm
Accept Reject Algorithm Implementation
In[224]:=
Test case for n=10,000
Test the above function, and make a plot of histogram overlaid on top of density of N(0,1)
In[225]:=
Out[232]=
In[233]:=
Out[233]=
The above is a plot showing the histogram for random numbers generated using the accept - reject method for N=10,000. The random numbers are very close the N[0,1] which indicates this method is working well. The larger N is, the more closely the random numbers histogram will approach N[0,1] probability density.
I have implemented a GUI based simulation as well for the above problem, please see the appendix below to run the simulation part.
Problem 3 simulation
In[234]:=
In[235]:=
Out[235]=
Created by Wolfram Mathematica 6.0 for Students - Personal Use Only (25 September 2007) |