##### Demonstration Script #0a
##### MATH322
#####
##### Brief introduction to (Monte Carlo) simulation
#####


### Introduction
set.seed(370)           # Ensure we agree on the numbers
x = runif(10)           # Draw 10 from a Unif(0,1) distribution
mean(x)                 # Calculate sample mean
sd(x)                   # Calculate sample sd
median(x)               # Calculate sample median




### What are the stats of Unif(0,1)?
set.seed(370)           # Ensure we agree on the numbers
x = runif(1e7)          # Draw 10 million from a Unif(0,1) distribution
mean(x)                 # Calculate sample mean (close to population mean)
sd(x)                   # Calculate sample sd (close to population sd)
median(x)               # Calculate sample median (close to population median)




### Another distribution: Student's t
set.seed(370)           # Ensure we agree on the numbers
x = rt(1e7, df=19)      # Draw 10 million from t(df=19) distribution
mean(x)                 # Calculate sample mean (close to population mean)
sd(x)                   # Calculate sample sd (close to population sd)
median(x)               # Calculate sample median (close to population median)




### Upshot:
# This method can be used to estimate the population statistics,
# but only if those values exist.
#
###





#####
### The Central Limit Theorem

# Standard Uniform distribution
mn=numeric()            # Set aside memory for the mn variable

for(i in 1:1e4) {       # Start for loop (length=10,000)
 x = runif(1000)        # Draw 1000 values from a Unif(0,1)
 mn[i] = mean(x)        # Calculate sample mean
}                       # End loop

hist(mn)                # Draw empirical pdf (aka histogram)

# Note that the histogram has the bell-shape to it. It is close to
# a Normal distribution, as the Central Limit Theorem specifies.




# Norm(1,100) distribution
mn=numeric()                 # Set aside memory for the mn variable

for(i in 1:1e4) {            # Start for loop (length=10,000)
 x = rnorm(1000, m=1,s=100)  # Draw 1000 values from a Norm(1,100) where 100 is sd, not var
 mn[i] = mean(x)             # Calculate sample mean
}                            # End loop

hist(mn)                     # Draw empirical pdf (aka histogram)

# Note that the histogram has the bell-shape to it. It is close to
# a Normal distribution, as the Central Limit Theorem specifies.




# X2(df=1) distribution
mn=numeric()                 # Set aside memory for the mn variable

for(i in 1:1e4) {            # Start for loop (length=10,000)
 x = rchisq(1000, df=1)      # Draw 1000 values from a X2(df=1)
 mn[i] = mean(x)             # Calculate sample mean
}                            # End loop

hist(mn)                     # Draw empirical pdf (aka histogram)

# Note that the histogram has the bell-shape to it. It is close to
# a Normal distribution, as the Central Limit Theorem specifies.




# Cauchy distribution
mn=numeric()                 # Set aside memory for the mn variable

for(i in 1:1e4) {            # Start for loop (length=10,000)
 x = rcauchy(1000)           # Draw 1000 values from a Cauchy dist
 mn[i] = mean(x)             # Calculate sample mean
}                            # End loop

hist(mn)                     # Draw empirical pdf (aka histogram)

# Note that the histogram does *not* have the bell-shape to it. This is because
# the Cauchy violates the requirement of the CLT, namely the finite variance 
# requirement.
#



###
# What does the Cauchy look like?

x=seq(-5,5,length=1000)
y=dcauchy(x)
plot(x,y)

# or, even better

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

# Note that it looks Normal, but it is not. There is too much "weight"
# in the tails.
#
# By the way, the Cauchy is also the t(df=1) distribution. So, as the df 
# increases from 1 to Inf, the t distribution changes from a Cauchy to
# a Normal. It is a very flexible distribution if you are investigating 
# kurtosis.


###
# What does the Unif(0,1) look like?

x=seq(-5,5,length=1000)
y=dunif(x)
plot(x,y, type="l")


###
# What does the Norm(0,1) look like?
x=seq(-5,5,length=1000)
y=dnorm(x)
plot(x,y, type="l")


###
# What does the X2(df=1) look like?
x=seq(0,5,length=1000)
y=dchisq(x, df=1)
plot(x,y, type="l")


###
# What does the F(1,4) look like?
x=seq(0,5,length=1000)
y=df(x, df1=1,df2=4)
plot(x,y, type="l")


###
# What does the Beta(1,3) look like?
x=seq(0,1,length=1000)
y=dbeta(x, shape1=1,shape2=3)
plot(x,y, type="l")


###
# What does the Beta(3,3) look like?
x=seq(0,1,length=1000)
y=dbeta(x, shape1=3,shape2=3)
plot(x,y, type="l")




### Upshot:
# There are several distributions we will be working with in 
# this course. These are the most helpful for certain cases.