Lab1: Introduction to R and Simulating Random Processes

PART I

Getting Started with R

Click on the R icon

If no R icon, go to the C drive and the Programs folder to find R.

To quit R:

> q()

What does it ask you? What does this mean?

Example 1: Let the games begin.

Generate a sample of 100 uniform random variables, Uniform(0,90).

For R HELP, for example to get help about the function runif:

> help(runif)

> runif(100, 0, 90)

What happens with the numbers?

It would be better to:

> temp <- runif(100, 0, 90)

What is in the object temp?

What is the length of the oject temp?

Make a informative plot of temp?

Example 2: The Binomial Distribution

For the binomial distribution, these functions are pbinom, qbinom, dbinom, and rbinom.

For help, use the help function on one of the four functions above.

Problem: Dr. Dribble has a probability of .8 of making free throws each time he shoots. What is the probability of him making at least 8 out of 10 free throws?

Assume his shots are independent of each other.

Let X be the number of free throws made. X is has a Binomial(10, 0.8) distribution. What does this look like?

Let's generate a large number of Binomial(n,p) random variables and look at the histogram. 

> set.seed(1)
> bindat <- rbinom(n=10000, size=10, prob=0.8)
> hist(bindat, breaks=seq(2, 10, 1), freq=F)
We want to find the probability of making at least 8 out of 10 free throws.

Let's calculate P(X <= 7) when X is has the Binomial(10, 0.8) distribution using the pbinom function.

> pbinom(7, 10, 0.8)

Now, to answer the question: P(X >= 8) = 1 - P(X <= 7) = 1-0.3222005 = 0.6777995

Example 3: The Exponential Distribution

For the exponential distribution, these functions are pexp, qexp, dexp, and rexp.

For help, use the help function on one of the four functions above.

Problem frow HW#1 Exam P: Let T be the number of days that elapse before a high-risk divier is involved in an accident. Assume T is exponentially distributed with rate λ=0.00713. What proportion of high-risk drivers are expected to be involved in an accident during the first 80 days of the calendar year?

Example 4: Let's make some plots in R!

1) In HW#2, there was a SDF that you were asked to plot. 3.4(b)
S(x) = 1/(1+x)^2, x>0.

Let's generate a sequence of x-values, evaluate the SDF function values, and plot. 

> x <- seq(0,10,by=.01)
#or x <- seq(0,10,length.out=1000)
> Sx <- 1/(1+x)^2
> plot(x, Sx, type="l")
> title("S(x) for 3-4(b)")

2) Let's try 3.4(a):

3) Let's try the Gompertz SDF and force of mortaility: B = 0.0003, c = 1.07 

#2x2 plotting region
> par(mfrow=c(2,2))
> x <- seq(0,100,by=.01)
> B <- 0.0003
> c <- 1.07
> Sx <- exp((B/log(c))*(1-c^x))
> plot(x, Sx, type="l")
> title("Gompertz SDF")
> mux <- B*c^x
> plot(x, mux, type="l")
> title("Gompertz Force of Mortality")
> plot(x, mux*Sx, type="l")

Example: Let's make some more plots in R!

Let's look at plotting different distributions

In order to run the function, you will need:

demo.gamma2.r

There are 2 ways to get these functions for your very own.

  1. Copy and paste the code into the R command window each function.
  2. Use the R source command which sources the code in the file:

> source("http://www.rohan.sdsu.edu/~babailey/stat575/demo.gamma2.r")

Now, to run the demo:

> demo.gamma2()

Other Probability Distributions in R

Let's Use Dr. Geyer's great webpage!
Work though this webpage and you can just copy and paste the R code in the Rweb box in R.
You do not need to use Rweb.

Here is the link

Intro to R and for more practice with R

Here is an extensive: Introduction to R

Short Intro to R:

If you want even more practice:
Try going through the following Intro. Lessons. The # sign is a comment, so you can copy and paste.

  • R vectors, lists and functions
  • R arithmetic and logical operators
  • R simple functions
  • R univariate random variable generators

    PART II

    Example 3: Markov Chain, Longrun Prabability and the Gamber's Ruin

    Let's look at my code in: game.r

    You can use this code for Homework 9!

    Using Matrices in R

    Return to Class Game Example:

    Q1: What is the probability that the gambler quits in six (or fewer bets)?
    Beginning with $20, the gambler has the following probability of being in
    an absorbing state at time 6: 0.616 + 0.274 = 0.899
    So, the gambler has an 89% chance of placing six or fewer bets!

    Let's make a P matrix and then calculate P^6.

    Q2: What is the expected number of bets he makes until he quits?

    A gambler who begins with $20 will bet an average of 3.85 times before either loosing all \ his money or doubling it!
    (A gambler who begins with $10 will bet an average of 2.54 times before he quits,
    and a gambler who begins with $30 will bet on average 3.31 times before he quits)

    Example 3: Gambler's Ruin Using R

    Here is a description of the simulation

    Let's check out the R code

    Example 4: Longrun Probability Using Java

    To simulate repeated game-playing, until one side or the other is ahead by a sufficient number,
    Let's check out longrun

    Example 5: Gambler's Ruin Using Java

    Let's check out Gambler's Ruin Web Page

    Here's another site: gambler

    PART III

    Example 6: Simulate a Markov Chain

    Let's look at the code in: mcsim.r

    What does this do?

    If you copy and past the code you should have a trans and run

    In order to run the run function you will need to:

    > run()

    Example 7: Cards and Poker

    We will need 3 functions to play poker. I have written these functions, so the help files are not in R, yet.

    They are not too complicated. The 3 functions are:

    cards.r
    rperm.r
    poker.r

    There are 2 ways to get these functions for your very own.

    1. Copy and paste the code into the R command window for each function.
    2. Use the R source command which sources the code in the files. You will need 3 source commands:

      > source("http://www.rohan.sdsu.edu/~babailey/stat550/cards.r")
      > source("http://www.rohan.sdsu.edu/~babailey/stat550/rperm.r")
      > source("http://www.rohan.sdsu.edu/~babailey/stat550/poker.r")

    For practice, generate 1 hand of 5 card poker.

    IMPORTANT If you want to reproduce your results, you will have to set the seed of the random number generator.
    There is an R function set.seed

    To use it,

    > set.seed(1)
    > poker(draw=5, hand=1)

    How can you draw the same hand again?

    Looping in R

    Let's play 10 poker hands:

    > rhands <- matrix(NA,10,5)
    > for (i in 1:10){ rhands[i,] <- poker(5) }

    What's rhands?

    Intro to R Part I: For more practice with R

    Here is a first introduction to R: ISU

    At the bottom of the page is how to get the "An Introduction to R"

    Part IV

    Generating Random Normal Variables - Old School

    Let's check out the Box-Muller Transformation!

    Here is a reference: Wolfram's page

    Example 8: Example 5.7 from text using Monte Carlo Integration

    An AC unit has an exponential lifetime distribution with mean 4 years.
    Use an effective annual rate of interest of 5%.
    Let Zbar denote the present value random variable for the unit of payment made at the time point of failure.
    Find E(Zbar). (the exact asnwer using integration by parts is 0.83670)

    Let's use Monte Carlo Integration to solve this problem!

    Let's take a look at a Lab from Dr. Hancock at Clark University

    What expectation (integral) do we need estimate?

    Here is my R code for this example. You can copy and paste the code into the R command window to run it.