######################
# 2014-04-24
# classical t-test
# example
######################


rm(list=ls())
par(mfrow=c(1,1))

source("scripts/loadano2.R")

#station<-"USW00014735"
station<-"USW00094728"
month<-"Jan"
buffer<-loadano2(station=station,month=month,start=1950,end=1980)
x<-buffer$ano


# Null hypothesis: the 1950-1980 temperature anomaly not different from zero


# test statistic
x0<-0
n<-length(x)
m<-mean(x)
vx<-var(x)
vm<-vx/n
# assuming independence from year to year in the monthly (or annual) temperature
# anomalies the variances of
d<-m-x0
t<-d/sqrt(vm)


print("sample size, sample mean, x0, difference mean-x0")
print(paste(n,round(m,4),round(x0,4),round(d,4)))
print("sample stddev, sample mean stddev, t-statistic, deg. freedom")

# if H0 is true then t follows a student t-distribution
# with n-1 degrees of freedom
df<-n-1
print(paste(round(sqrt(vx),4),round(sqrt(vm),4),round(t,4),df))

# p-value from the cumulative density distribution

pvalue<-pt(t,df=df)
print(paste("Prob(t<=",round(t),"given Ho was true)=",round(pvalue,8),sep=''))

# built-in student t-distribution
# dt
thelp<-seq(-6,6,0.1)
pdf<-dt(thelp,df=df)
cdf<-pt(thelp,df=df)
plot(thelp,pdf,typ='l',lty=2,col="black",
     main='Student T-test',xlab='t value',ylab='PDF (and CDF)',
     ylim=c(0,1)
     )
lines(thelp,cdf,col=1,lty=1)
points(c(t,t),c(0,0),col="green",pch=19)
lines(c(t,t),c(0,pvalue),col='green')