### Statistics 6401, Handout 1

### This R script simulates flipping a fair coin 1000 times.
### Three examples of using simulation are given using various 
### R commands.  Finally, the CLT is demonstrated.

# Example 1. Simulate flipping a fair coin n times.

n = 1000
p = 0.50

h = 0 	# counter for heads

for (i in 1:n){	
	x = runif(1)
	if (x < p) h = h + 1
}
h	# total number of heads simulated

ph = h/n
ph	# simulated estimate of the probability of heads

# Example 2. Simulate flipping a fair coin n times.
# Simulate using vectors.

x = runif(n)	# creates a vector x of random uniform values

length(x)

x[1:10] # prints the first 10 values of the vector x

C = x[x < p]	# puts all values of x that are less than p into a vector C
			
h = length(C)	# gives the number of heads
h

ph = h/n
ph

# Example 3. Recall that flipping a fair coin n times is a binomial experiment
 
h = rbinom(1,n,p)
h

ph = h/n
ph

# Example 4. Now simulate flipping a fair coin n times, but do it 100,000 
# times and then plot plot the sampling distribution of the estimated 
# probability of heads.  What does the CLT say about the distribution?  
# Ans: approximately N( p, sqrt((p*(1-p))/n) )

Reps = 100000

h = rbinom(Reps,n,p) # this produces a vector of length Reps

ph = h/n	# this uses the vector division elementwise

ph.min = min(ph)		# these commands find the range of values of ph
ph.max = max(ph)

hist(ph, prob=T, xlim=c(ph.min,ph.max), main="Sampling distribution of the probability of heads")

par(new=T) # resets the plot so another plot can be made of the same graph sheet

x = seq(ph.min, ph.max, 0.001)

p.mean = p
p.sd = sqrt((p*(1-p))/n)

y = dnorm(x, p.mean, p.sd)

plot(x, y, type="l")

####################################################################

### This is an R program that simulates a system of 
### n components connected in series.  The probability that
### a component fails is p.  The probability that
### the system fails is approximated using simulation.

p = 0.25 
n = 4

Reps = 10000

Count2 = 0	# Counter for number of failing systems in the overall simulation
for (j in 1:Reps){
	Count1 = 0	# Counter for number of failures in a simulated system
	for (i in 1:n){
		ci = runif(1)
		if(ci < p) Count1 = Count1 + 1
	}
	if (Count1 > 0) Count2 = Count2 + 1
}

psf = Count2/Reps
psf

psf.truth = 1-(1-p)**n
psf.truth

### Alternative solution, using R

c = rbinom(Reps,n,p)
Count = c[c>0]
psf = length(Count)/Reps
psf

psf.truth

### Exercise.  Write and R program that simulates a system of 
### n components connected in parallel. The probability that
### a component fails is p.  

####################################################################

### This script plots the normal density four times with different 
### parameter values.

par(mfrow=c(2,2)) # makes a grid of plots

# plot 1

mu = 0	# standard normal
sigma = 1

x.min1 = mu - 3*sigma
x.max1 = mu + 3*sigma

x1 = seq(x.min1,x.max1,0.01)
y1 = dnorm(x1,mu,sigma)

plot(x1,y1,type="l", main="Normal(0,1)")

# plot 2

mu = 1	# standard normal
sigma = 1

x.min2 = mu - 3*sigma
x.max2 = mu + 3*sigma

x2 = seq(x.min2,x.max2,0.01)
y2 = dnorm(x2,mu,sigma)

plot(x2,y2,type="l", main="Normal(1,1)")

# plot 3

mu = 1	# standard normal
sigma = 2

x.min3 = mu - 3*sigma
x.max3 = mu + 3*sigma

x3 = seq(x.min3,x.max3,0.01)
y3 = dnorm(x3,mu,sigma)

plot(x3,y3,type="l", main="Normal(1,2)")

# plot 4, plot all densities on the same plot

x.min = min(x.min1,x.min2,x.min3)
x.max = max(x.max2,x.max2,x.max3)

y.max = max(y1,y2,y3)

plot(x1,y1,type="l",lty=1,xlim=c(x.min,x.max),ylim=c(0,y.max),main="Normal Densities",xlab="x",ylab="y")
par(new=T)
plot(x2,y2,type="l",lty=2,xlim=c(x.min,x.max),ylim=c(0,y.max),xlab="x",ylab="y")
par(new=T)
plot(x3,y3,type="l",lty=3,xlim=c(x.min,x.max),ylim=c(0,y.max),xlab="x",ylab="y")

### Exercise.  Make plots for the exponential, gamma, t and F distributions.