#Simple session with R -- Lecture 1

#checking directories

dir()

dir("/Users/mats-ksp/Dropbox/Public/Webpage/Teaching/DataMining/Data")

setwd("/Users/mats-ksp/Dropbox/Public/Webpage/Teaching/DataMining/Data")

dir()

#Getting data in a vector
x=scan("Table2_1.txt")

#The following would work too

xx=scan("/Users/mats-ksp/Dropbox/Public/Webpage/Teaching/DataMining/Data/Table2_1.txt")

n=length(x)
mean(x)
sd(x)


#Bootstrapping means
help(sample)
B=100
Bmean=vector('numeric',100)
for(i in 1:B)
{
Bmean[i]=mean(sample(x,n,rep=T))
}
sd(Bmean)

sd(x)/sqrt(n)

#Bootstrapping variances
B=1000
Bvar=vector('numeric',B)
Bsd=Bvar
for(i in 1:B)
{
Bvar[i]=var(sample(x,n,rep=T))
Bsd[i]=sqrt(Bvar[i])
}
sd(Bvar)

sd(Bsd)

#Data 
n=100
x=rnorm(n,5,10)

#Bootstrap
B=10000
Bvar=vector('numeric',B)
for(i in 1:B){
Bvar[i]=var(sample(x,n,rep=T))}
mean(Bvar)
sd(Bvar)
hist(Bvar,nclass=10)

#Monte Carlo
N=15000
MCvar=vector('numeric',N)
for(i in 1:N){MCx=rnorm(n,5,10)
MCvar[i]=var(MCx)}
mean(MCvar)
sd(MCvar)
X11() #graphical window in Unix
      #windows() in Windows
      #quartz() in Mac
hist(MCvar,nclass=10)


#Non-linear regression

X=runif(50,0.5,8)
e=rnorm(50,0,0.35)
Y=sin(X)+e

plot(X,Y)


#EXTRA CODE NOT USED IN 2014 LECTURE!
#BEGINING
#Parametric bootstrap -- simple example

#Data

x=scan("Table2_1.txt")
n=length(x)
mean(x)
sd(x)

#Bootstrapping variances

B=1000
Bvar=vector('numeric',B)
Bsd=Bvar
for(i in 1:B)
{
Bvar[i]=var(sample(x,n,rep=T))
Bsd[i]=sqrt(Bvar[i])
}
sd(Bvar)

sd(Bsd)

hist(Bvar,nclass=10)

#Monte Carlo study of variances

N=15000
MCvar=vector('numeric',N)
for(i in 1:N){
MCx=rnorm(n,50,0.1)
MCvar[i]=var(MCx)
}
mean(MCvar)
sd(MCvar)
X11() #graphical window in Unix
      #windows() in Windows
      #quartz() in Mac
hist(MCvar,nclass=10)


#Parametric bootstrap

#Fit the model -- normal distribution

#Data

x=scan("Table2_1.txt")
n=length(x)
mu=mean(x)
sigma=sd(x)

#Simulate from the fit

PB=1000
PBvar=vector('numeric',PB)
PBsd=PBvar
for(i in 1:PB)
{
PBx=rnorm(n,mu,sigma)
PBvar[i]=var(PBx)
}

mean(PBvar)
sd(PBvar)
X11() #graphical window in Unix
      #windows() in Windows
      #quartz() in Mac
hist(PBvar,nclass=10)

#END