# -------------------------------------------------------------------
# Lecture 3: Normal distribution, estimators, sampling distributions of estimators, tests of hypotheses

# Required libraries: --
  rm(list=ls())
  # source("http://klein.uk/R/myfunctions.R")
  ls()
# -------------------------------------------------------------------


# --- Standard normal distribution ---------------------------------

# Normal distribution: 68-95-99.7 or six sigma rule
  plot(dnorm(seq(-5,5,.1)) ~ seq(-5,5,.1), type="l")
  abline(h=0, v=c(-3,-2,-1,0))

# What proportion of observations are smaller than 0.83? (p.6)
  pnorm(0.83)

# What proportion of observations are greater than -2.15? (p.7)
  1 - pnorm(-2.15)

# Inverse of SND: F^{-1}(0.03) = ? (p.8)
  qnorm(0.03)

# Inventory example. (p.9ff)
  # What is the probability of a stockout?
  1 - pnorm(20, mean=15, sd=6)
  # If the prob of stockout is to be no more than 5%, what should the reorder point be?
  qnorm(1-0.05, mean=15, sd=6)




# --- Asymptotic properties of estimators (pages 23-25) ----------------

  # Simulate probability distribution of sample mean.
  d = function(n) density( sapply(1:10000, function(x) mean(rnorm(n, mean=100, sd=50))) )

  # Plot probability density of sample mean (sample size 1).
  plot( d(1), ylim=c(0,.08), main="Distribution of sample mean")
  abline(h=0)

  # Sample size 4, 25, and 100.
  sapply(c(4,25,100), function(x) lines( d(x) ) )




# --- Simulation: Sample variance, unbiased estimators (pages 27-29) ----

# --- PART A: Biased vs unbiased estimator for population variance.

# What's special about n-1 in the equation for the sample standard deviation?
# What would happen if we used n instead?

# Unbiased sample variance with denominator n-1
  var(0:1)

# Biased sample variance with denominator n
  newvar = function(x) 1/length(x) * sum( (x-mean(x))^2 )
  newvar(0:1)

# 1,000 Simulations (sample size 10) and mean of biased vs unbiased sample variance
  s2 = sapply(1:1000, function(x){
      sample10 = rnorm(10,mean=0,sd=1)
      c( var(sample10), newvar(sample10) )
    }
  )
  s2.unbiased = s2[1,]; mean(s2.unbiased)
  s2.biased = s2[2,]; mean(s2.biased)


# --- PART B: Consistency of the biased variance estimator.

# If the sample size is increased it ceases to matter whether we use n or n-1 in the denominator.
# Simulate distribution of biased sample variance for different sample sizes
  d = function(n) density( sapply(1:10000, function(x) newvar(rnorm(n, mean=100, sd=50))) )

  # Sample size 10
  plot(d(10), ylim=c(0,.005)); abline(h=0, v=50^2)

  # Sample sizes 20, 100, and 1000
  sapply(c(20,100,1000), function(x) lines(d(x)) )

# In summary, using n gives a biased estimate of the true variance. The smaller the sample size, the greater this discrepancy between the unbiased and biased estimator.




# --- Simulation: Power of a test (pages 67-69) -------------------------

# Define hypothesis test function.
  h.test.A = function(n, mu){
    sapply(mu, function(x) abs( mean( rnorm(n, mean=x, sd=1) ) ) > qnorm(0.975)/sqrt(n) )
  }

  h.test.B = function(n, mu){
    sapply(n, function(x) abs( mean( rnorm(x, mean=mu, sd=1) ) ) > qnorm(0.975)/sqrt(x) )
  }


# --- PART A: Power of a test and evidence against H0 -------------------

# Set values of sample size and population mean.
  n = 10
  mu = c(0, .05, .1, .2, 1, 2)
# Run simulation for sample size n=10 and population means of 0, .05, .1, .2, 1, and 2.
  data = sapply(1:10000, function(x) h.test(n=n, mu=mu))
# calculate percentage of rejections when null is not true (= power of test).
  rejections = sapply(1:6, function(x) sum(data[x,])/10000)
# To see the test's power, graph the prob of rejecting H0 against the evidence.
  plot(rejections ~ mu, ylab="Prob of rejecting H0", xlab="Evidence against H0", main="Power of test")
  lines(lowess(rejections ~ mu, f=0.5))
  abline(h=0.05)  # horizontal line for size of the test, i.e., Prob(rejecting H0 | H0 true)


# --- PART B: Power of a test and sample size ---------------------------

# Set values of sample size and population mean.
  n = c(10, 100, 1000)
  mu = 0.2
# Run simulation for population mean mu=0.2 and sample size of 10, 100 and 1000. 
  data = sapply(1:10000, function(x) h.test.B(n=n, mu=mu))
# calculate percentage of rejections when null is not true (= power of test).
  rejections = sapply(1:3, function(x) sum(data[x,])/10000)
# To see the test's power, graph the prob of rejecting H0 against the evidence.
  plot(rejections ~ n, ylab="Prob of rejecting H0", xlab="Sample size", main="Power of test")
  lines(lowess(rejections ~ n, f=0.5))




# -------------------------------------------------------------------
# --- End of Session ------------------------------------------------

  q("no")