Fall Semester 2000
This is an exercise to produce non-uniform variates from a standard uniform random number generator.
Things looks bad for Alice's CryptoPhone company. Her chief analyst, Bob, has just discovered that the microfleem controllers in her customers' secure-cell phone units are defective. In fact, the microfleems are starting to fail purely at random, with an average time to failure of 3 months. The distribution of time-to-failure values for a pure random process is given by an exponential,
p(t) = (1/tau)*exp(-t/tau) [for t>=0, otherwise 0]
where p(t) dt is the probability that a microfleem will fail
at some time between t and t +dt; the mean failure time,
tau, is 3 months.
In this exercise, write a C++ code called rndexp.cc that produces a list of simulated microfleem failure times relative to the present time. Use time units of months. When your code is compiled and linked, it should behave as follows:
darkstar> rndexp 1000 3.0 314159265
3.302108784
4.753503821
...
so that your code prints a column of 1000 exponentially distributed,
random numbers with a mean value of 3.0, and with 314159265 as a seed
for the "random" sequence.
A good starting place for this exercise is your "random.cc" code from a04.html, exercise 2. Feel free to use this solution code in directory ~p5720/solutions/a04. There are some differences in details between random.cc and the code you are to write for this exercise: Your new code should take three command line arguments instead of two. The (new) second argument, the mean of the exponential distribution, should be of type double (thus the "ascii-to-double" function atof() may come in handy, as did the "ascii-to-int" function atoi()). Next, you should get rid of "myrandom()" and "myrandomseed(seed)" (both the function definitions below main(), and the function prototypes above main(), and instead generate random numbers with the Mersenne twister in function "drand()" (this returns a uniform random number in double precision). To seed drand(), use "sdrand(seed)", where seed is any positive integer. Although source code is available in ~p5720/src/random.cc (sorry about the name conflict with the code from Assgn. 4), you should not use it in your code. Instead, use the class C++ library -- it contains drand() and sdrand(). To access the p5720 library, first #include the "p5720.h" header file (located in ~p5720/include) in your source code like this:
#include "p5720.h"This header has function prototypes for all functions in the p5720 class library, including drand() and sdrand(). The compiler can then link to the library itself using the "-lp5720" option. In this instance, neither the header or the library archive file are located in standard places. (They are tucked away in directories within ~p5720/.) Thus we need to feed the compiler special directives using "-I" (for #iinclude) and "-L" (for Library) options to indicate where it should look for the header and library. For example,
g++ -I/u/course/p5720/includetells the compiler to find the header "p5720.h" in directory ~p5720/include.
This compiler option is cumbersome to type in, so you can avoid having to do so with the "make" utility. "Make" uses "makefiles" to give it instructions on how to build executable programs (among other things). It can be easily set up to handle longwinded compiler options like -I. (It also is great for organizing compilations that involve a large number of separate source codes.)
A makefile for this assignment that takes care of the details of the -I, -L and -l options, is located in ~/p5720/makefiles/a07.mk. You can use it by typing
make -f ~p5720/makefiles/a07.mk rndexp
This will compile your rndexp.cc code and create an executable program
"rndexp". The -f option tells which makefile to use.
If you prefer not to have to type in the "-f ~p5720/.../a07.mk", you can
simply copy a07.mk to your working directory, calling it "Makefile".
Then just type
make rndexp
(When no -f option is given, "make" by default looks for a file called
"Makefile" (or "makefile") to figure out what to do for the specified
target (here, "rndexp").
As for designing the algorithm to generate exponentially distributed random numbers, recall that you can use a uniform random number generator pumped through some function (specifically the inverse of the cumulative distribution function). To do this, follow the presciption given in lecture which was in turn based on the write up in Numerical Recipes (Section 7.2). (You can obtain NR text on-line at this URL -- click on the Books On-Line icon.) To get a unit-mean exponential distribution, start with the probability density function p(x) and its cumulant P(x),
p(x) = exp(-x) [for x>=0, 0 otherwise]
P(x) = Integral[0 to x] dz*p(z) = 1-exp(-x).
Then invert the cumulative distribution P(x), to find a function invP(y)
such that invP(P(x))=x. In this case you can prove that InvP(y)= -log(1-y)
[just substitute y=1-exp(-x)]. Finally, realize that 1-y is a number between
0 and 1 and that you can get a exponential random x-value simply by setting
"y" in invP(y) to "u", a uniform random number on [0,1) as produced by
drand(), for example. Equivalently, you can use (1-y)->u to get an even
simpler expression:
x_random= - log(drand());
Here you only have to figure out how to do the same thing
with the probability density function
p(x) = 1/tau * exp(-x/tau);
with tau, the mean of the distribution, set to 3.0.
Hints and Comments:
setenv ARCHTYPE sparc # if you are on a Sun (dark/blackstar or fusion) setenv ARCHTYPE rs6000 #if you are on an IBM (gp##)
Turn in only your C code, rndexp.c.
This exercise is to write a very short Bourne (sh) shell script. The problem at hand is a continuation of exercise 1.
Through a clever multi-unit distributed cross-talk scheme, Alice finds a way to keep all of her cell phones operating for a while, even after some of the microfleems fail. Unfortunately, when 60 percent of the microfleems fail, the scheme will also break down. For how many months after the present time will Alice's scheme work?
To answer this question, write a Bourne shell script called fleemtime which uses your rndexp program from exercise 1 to produce output like this:
blackstar> fleemtime
Months to 60 percent fleem failure: XXX
You can start with this script:
#!/bin/sh # a starting point for fleemtime: answer=`rndexp 1000 3.0 314159 | head -1` echo "The variable answer is $answer"You can put shell commands in here, like "echo", or you can set variables like "answer". Note that the single backquotes (`) mean run the quoted command and turn the result into a string. (Check out the example shell scripts in ~p5720/examples, maybe looking at sh-script-simple first.)
Hint 1: The utilities head, tail and sort may be useful.
Hint 2: To run your script like an executable, be sure to
give at yourself permission to execute, using the Un*x "chmod" utility:
chmod 700 fleemtime
("700" implies read-write-execute for you, no privileges for anyone else.)
Submit your script fleemtime.
This exercise is a first example of Monte Carlo integration. (We will do another later on.)
The active fault regions in the Wasatch Front experience an earthquake of magnitude greater than 7 on the Richter scale every 350 years, on average. The goal of this exercise is to estimate the probability of an earthquake occurring in one of the active regions, the SLC fault zone. It is known that the last quake of magnitude greater than 7 on the SLC fault zone occurred roughly 2600 years ago, and that the average time between quakes is around 1600 years. By assuming a specific model for the distribution of time between quakes, you can calculate the probability that there will be an earthquake of magnitude 7 or larger within, say, the next 50 years.
In this exercise, assume that the time between quakes, T, is distributed lognormally, that is, the logarithm of the variate T is Gaussian. Lognormal distributions are characteristic of some random quantities which cannot have negative values. A famous example is the price of stocks (this is the basis for the Black-Scholes equation, which led to the 1997 Nobel prize in Economics). Lognormal models for earthquakes (or even quakes on neutron stars) are popular among researchers because they are broadly peaked yet give vanishingly small values of earthquake probabilities immediately following a major quake. This fits in well with the idea that a major quake releases geological stresses, and a second major quake can follow only after these stress can build up again.
Here, write a C++ code, quake.c++, that simulates the time between quakes so that you can determine by multiple trials the likelihood that a magnitude 7+ quake will occur within the next 50 years. For definiteness, use parameters for an SLC earthquake as given in an article by McCalpin and Nishenko [1996 (Feb), Journal of Geophysical Research]:
S = sqrt(log(sd*sd/(mu*mu)+1.0));
M = log(mu)-0.5*S*S;
x = exp(S*gaussrand()+M);
As in Exercise 1 above, use the makefile
"a07.mk" to compile and link your quake.cc code. The makefile
will automatically connect you with the p5720 header and library
where "sdrand()" and "gaussrand()" are declared/defined.
Use target name "quake" when calling the "make" utility (that
is, just substitute "rndexp" with "quake" in the "make" commands
given in Ex. 1).
Ultimately, your code should run exactly like
blackstar> quake 10000 987654321
percent chance of quake in 50 years: XXX
where 10000 is the number of trials and 987654321 is a seed
for the random number generator.
(Again, a good starting point for your code might be the
random.cc program from
Assignment 4, Exercise 2. However, be sure to
call "sdrand()" before "gaussrand()", and do not
print out anything from within your loop over trials!)
Hint: When you test whether a random time-to-quake is between now and 50 years from now, don't forget that this range corresponds to 2600 to 2650 years since the last earthquake. The algorithm described here simply produces lognormal samples, even though most of the samples give quakes that occur before today (t = 2600 years). The vast majority of trials which yield "early" quakes should be completely ignored. Since your answer may depend on a relatively small subset of the total number of trials, you may wish to vary the [perhaps large] number of trials to roughly estimate the minimum value of "n" needed to get a stable result. As a check to see if your algorithm is correct, consider that the chance of an earthquake from now until t = many, many (>> 50) years from now is 100%.
Submit your C++ code quake.cc.