#### MEETING 5
####
#### NORMAL DISTRIBUTION
#### CONFIDENCE SETS
####


#### DATA & COMMON FUNCTIONS
#### http://people.ucsc.edu/~mwagers/ling280/data/chamorro.Rdata

source("http://people.ucsc.edu/~mwagers/ling280/notes/ling280.R")
load(url("http://people.ucsc.edu/~mwagers/ling280/data/chamorro.Rdata"))

#### Two data frames - judgment times in classifying good/bad extractions
#### poss$RT / RTs to possessor extractions
#### arg$RT  / RTs to verb argument extractions


###
### Normal distribution
###

my.norm.pdf <- function(...) {
  # 
  #
  # ARGS:
  #
  # RETURNS:

  a <- sigma * sqrt(2*pi)
  b <- (-1/2) * (( x - mu ) / sigma)^2

  fx <- a^-1 * exp(b)

  return(...)

}

x <- seq(-5, 5, by=0.01)

plot(x, my.norm.pdf(x, 0, 1), type='l', lwd=3, col="red")
hist(rnorm(1e6), freq=FALSE, breaks=1e2, add=1, density=15)

# 
# CDF (Phi)
#

my.norm.cdf <- function(x = 0, mu = 0, sigma = 1){
  # 
  #
  # ARGS:
  #
  # RETURNS:

  Fx <- integrate(my.norm.pdf, -Inf, x, mu, sigma)

  return(Fx$value)

}

my.norm.cdf()
my.norm.cdf(x = 5, mu = 5, sigma = 1)

my.norm.cdf(x = 6, mu = 5, sigma = 1)
my.norm.cdf(x = 6, mu = 5, sigma = 2)
my.norm.cdf(x = 7, mu = 5, sigma = 2)

my.norm.cdf(x = (7 - 5)/2 )
my.norm.cdf(x = (6 - 5)/1 )


# 
my.norm.cdf(1) - my.norm.cdf(-1)


# CDF 1
mu<-100
sigma<-10
n.obs<-1e4

rnorm(n.obs, mean = mu, sd = sigma) -> norm.1
ecdf(norm.1) -> norm.1.ecdf;

seq(-5, 5) -> x.sd
min.x <- mu - 5 * sigma 
max.x <- mu + 5 * sigma

plot(x.sd, 
  norm.1.ecdf(seq(min.x, max.x, by=sigma)),
  xlab="SD",ylab="P(X<x)", type='b')

lines(my.norm.cdf(100,10))

####
#### CONFIDENCE INTERVALS(/SETS)
####

#
rnorm(10000, 0, 1) -> norm.standard
quantile(norm.standard, c(0.005,0.025,0.05,0.165,0.5,0.835,0.95,0.975,0.995)) -> norm.quants
print(round(norm.quants,digits=2))

#
rnorm( 100, 50, 25 ) -> norm.4
se.n4 <- sd(norm.4)/sqrt(length(norm.4))
ci.95 <- c( mean(norm.4) - 2 * se.n4 ,
			mean(norm.4) + 2 * se.n4 )


######################
### T DISTRIBUTION  ##
######################

#
qt(0.025, df = 1:50)

rnorm( 20, 50, 25 ) -> norm.5
#
qt(c(0.025,0.975), df = 20 - 1) -> t.crit
#
qnorm(c(0.025, 0.975))

### CONFIDENCE INTERVAL
se.n5 <- sd(norm.5)/sqrt(length(norm.5))
ci.95 <- mean(norm.5) + t.crit * se.n5			

####
t.test(norm.5)

### 
### Are possessor extractions harder to process than argument extractions?
###

with(poss, tapply(RT, Subj, mean)) -> poss.subjectRTs
with(arg, tapply(RT, Subj, mean)) -> arg.subjectRTs