Lab1: Introduction to R

Lab1: Introduction to R and Simulating Random Processes

PART I

Getting Started with R

Click on the R icon

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: 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(c):

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")

NOW FOR MARKDOWN!

Intro to Markdown:

Let's look at Example/Problem from above Dr. Dribble:
Suppose the HW problem wants you calculate the probabilty AND also include a histogram.

Here is the code that we will need to put in a Markdown document: ExampleLab1.Rmd - Please see Canvas!
We will see if we can complete the Extra Problem (not to be handed in) for HW#4.

Here are some Rmarkdown resources
(Send me more and I will gladly add to this!)

NOW FOR SWEAVE!

Intro to Sweave:

Let's look at Example/Problem from above Dr. Dribble:
Suppose the HW problem wants you calculate the probabilty AND also include a histogram.

Here is the code that we will need to put in a Sweave document: drdribble

We will need to select R Sweave under the File, then New File options.

This automatically adds the lines:
\SweaveOpts{concordance=TRUE}

Here are some knitr resources
(Send me more and I will gladly add to this!)

EXTRA EXAMPLES - BELOW!

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.

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

> source("https://edoras.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

Here is a SHORTER Introduction to R Tutorial with Videos

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 1: Matrices in R and Markov Chains

Let's look at my code in: HW6 Problem 1 and 2

First let's just get the TPM matrix, and call it P

Let's make the P^2 and P^3. There are many ways to do this!

You can extract elements of the matrix, if needed! Can you extract the (3,3) element? The 3rd row?, etc.
You can also perform many other matrix operations (see Examples)

Now, back to HW6 and more! Let's make the Q matrix, the T matrix, and the QT matrix. (Hint: These might be helpful for practice to check your answers.)

We also know from lecture that the from QT matrix will can calculate the probability of absorption.

Example 2: Markov Chains R Package

Let's look at the: Crash Course Introduction to markovchain R package

Code in the above document is also available here.

First let's just get the TPM matrix, and call it P

Let's make the P^2 and P^3. There are many ways to do this!

You can extract elements of the matrix, if needed! Can you extract the (3,3) element? The 3rd row?, etc.
You can also perform many other matrix operations (see Examples)

Now, back to HW6 and more! Let's make the Q matrix, the T matrix, and the QT matrix. (Hint: These might be helpful for practice to check your answers.)

We also know from lecture that the from QT matrix will can calculate the probability of absorption.

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

Class Game Example:

A gambler bets $10 a game on a sequence of games. He starts with $20 and places bets one at a time until he either doubles ($40) his money or loses ($0) it all. Suppose he has a probability of 0.6 of losing and a probability of 0.4 of winning each time.

The amount of money he has after each game is a Markov chain with a one-step transition probability matrix. There are 5 States where $0 is State 1, $10 is State 2, $20 is State 3, $30 is State 4 and $40 is State 5. This Markov Chain has two absorbing states, $0 is State 1 and $40 is State 5.

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

(You can use this code for Homework!)

Using Matrices in R

Return to Class Game Example:

Q1: What is the probability that the gambler quits in six (or fewer bets)?

Ans: Let's make the P matrix and then calculate P^6.
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!

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

Ans: Using the Q matrix.
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)

Q3: What is the probability of absorption? What does this mean?

Example 3: The Gambler's Ruin

Here is the WOLFRAM Demonstration Project simulation

The description is below the simulation.

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.

Copy and paste the code into the R command window for each function.
Use the R source command which sources the code in the files. You will need 3 source commands:
> source("https://edoras.sdsu.edu/~babailey/stat550/cards.r")
> source("https://edoras.sdsu.edu/~babailey/stat550/rperm.r")
> source("https://edoras.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

Many tutorials available!

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

Here is a reference for Monte Carlo Methods: Monte Carlo Integration 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.

Example 9: Let's Bootstrap!

We're going to follow Lab: Nonparametric Bootstrapping Example 1