# Rozdział 6

#6.1

#Realizacja procedury bootstrap

data(cars)
attach(cars)
N=1000
z=dist/speed
n=length(z); wynik=matrix(0,N,1)
for (i in 1:N)
	{wynik[i,1]=mean(sample(z,n,replace=TRUE))}
print('Ocena wartości przeciętnej ilorazu') 
mean(wynik)
print('ocena wariancji estymatora') 
var(wynik[,1])

#6.2

#Konstrukcja przedziału ufności metodą bootstrap

x<-rbeta(30,2,2)
y<-rgamma(30,4)
z=cbind(x,y)
B=c()
for (i in 1:1000)
{
	n=sample(1:30,replace=TRUE)
	B[i]=sum(x[n])/sum(y[n])
}
### przedział ufności
quantile(B,c(0.025,0.975))


#6.3

#Konstrukcja przedziału ufności metodą bootstrap z wykorzystaniem pakietu boot

library(boot)
ratio <- function(d, w)sum(d$x * w)/sum(d$u * w)
boot_city=boot(city, ratio, R=1000, stype="w")
boot.ci(boot_city,type=c("norm","perc"))


#6.4

#Ocena wariancji estymatora metodą zależnych grup losowych

data(chickwts)
hist(chickwts$weight) 
x=sample(chickwts$weight,70)
z=c(0,0)
N=10
n=70/N
for (i in 1:N)
{
	z[i]=mean(x[((i-1)*n+1):(i*n)])
}

mean(z)

var1=var(z)/N
z0=mean(chickwts$weight)
var2=sum((z-z0)^2)/(N*(N-1))
print('Odchylenia standardowe - oceny')
sqrt(var1)
sqrt(var2)


#6.5

#Ocena wariancji estymatora metodą niezależnych grup losowych

data(chickwts)
hist(chickwts$weight)
x=sample(chickwts$weight,10,replace=N)
hist(x)
z=c(0,0)
N=100
for (i in 1:N)
{
	x=sample(chickwts$weight,10,replace=N)
	z[i]=mean(x)
}
mean(z)
var1=var(z)/N
z0=mean(chickwts$weight)
var2=sum((z-z0)^2)/(N*(N-1))
print('Odchylenia standardowe - oceny')
sqrt(var1)
sqrt(var2)


#6.6

#Porównanie oszacowań metodami bootstrap, DRG I IRG

x=rnorm(1000,100,5)
n=50; G=5
### metoda bootstrap
B=10000
sx=sample(x,n)
zB=c()
for (i in 1:B)
	{zB=c(zB,mean(sample(sx,n,replace=TRUE)))}
var(zB)
### metoda DRG
sx=sample(x,n)
zDRG=c()
for (g in 1:G)
	{zDRG=c(zDRG,mean(sx[(1+(g-1)*10):(g*10)]))}
var(zDRG)/G
### metoda IRG
zIRG=c()
for (g in 1:G)
	{zIRG=c(zIRG,mean(sample(x,10)))}
var(zIRG)/G


#6.7

#Symulacyjne porównanie oszacowań metodami bootstrap, DRG I IRG

wB=c(); wDRG=c(); wIRG=c()
for (k in 1:100)
{
	#= = == = = = =
	#Tu kompletny kod 6.6
	#= = == = = = ==
wB=c(wB,var(zB))
wDRG=c(wDRG,var(zDRG)/G)
wIRG=c(wIRG,var(zIRG)/G)
}
sd(wB);sd(wDRG);sd(wIRG)
par(mfrow=c(3,1))
hist(wB,xlim=c(0,2))
hist(wDRG,xlim=c(0,2))
hist(wIRG,xlim=c(0,2))


#6.8

#Procedura testu permutacyjnego - porównanie dwóch populacji

x1<-c(0,4); x2<-c(3,5,7)
x<-c(x1,x2); N=10000
T0=mean(x1)-mean(x2)
T=c()
for (i in 1:N)
{
	xs=sample(x)
	T[i]=mean(xs[1:2])-mean(xs[3:5])
}
hist(T)


# 6.9

#Kod procedury testu permutacyjnego z wykorzystaniem pakietu coin

library(coin)
x=c(0,4,3,5,7)
y=factor(c(1,1,2,2,2))
oneway_test(x~y)


#6.10

#Ocena prawdopodobieństwa zdarzenia rzadkiego

N=10^3
y=rexp(N)+4.8
waga=dnorm(y)/dexp(y-4.8)
sum(waga)/N
### Graficzna prezentacja wyników
plot(cumsum(waga)/1:N, type=”l”)
abline(a=pnorm(-4.8),b=0, col=”red”)


#6.11

#Generowanie wartości losowych z rozkładu normalnego z wykorzystaniem algorytmu Metropolisa-Hastingsa

n=1000
x=c()
x[1]=0
a=1
for (i in 1:n)
{
s=2*a*runif(1)-a
x[i+1]=x[i]+s
while (dnorm(x[i+1])/dnorm(x[i])<runif(1)) {s=2*a*runif(1)-a; x[i+1]=x[i]+s}
}
par(mfrow=c(1,2))
plot(x,type="l")
hist(x)


#6.12

#Generowanie wartości losowych z rozkładu chi-kwadrat z wykorzystaniem algorytmu Metropolisa-Hastingsa

n=1000
x=c()
x[1]=2
a=1
for (i in 1:n)
{
s=2*a*runif(1)-a
x[i+1]=x[i]+s
while (dchisq(x[i+1],3)/dchisq(x[i],3)<runif(1)) {s=2*a*runif(1)-a; x[i+1]=x[i]+s}
}
par(mfrow=c(1,2))
plot(x,type="l")
hist(x)


#6.13

#Generowanie wartości z rozkładu normalnego z wykorzystaniem biblioteki MCMCpack

library(MCMCpack)
wynik <- MCMCmetrop1R(dnorm, theta.init=5.3, mean=0,sd=1,thin=1,mcmc=10000,burnin=500,tune=2.0,verbose=5000,logfun=FALSE)
plot(wynik)


#6.14

#Generowanie wartości z rozkładu dwuwymiarowego na podstawie rozkładów warunko-wych – algorytm Gibbsa


t00=0.998; t01=0.121
y00=0.956; y01=0.080
ile=10000
y=c(); y[1]=0
t=c()
for (i in 2:ile)
{
	if (y[i-1]==1) t[i-1]=rbinom(1,1,1-t01) 
		else t[i-1]=rbinom(1,1,1-t00)
	if (t[i-1] == 1) y[i]=rbinom(1,1,1-y01)
		else y[i]=rbinom(1,1,1-y00)
}
table(y[1000:10000],t[1000:10000])

#6.15

#Zastosowanie algorytmu Gibbsa

N=1000
lambda=matrix(0,N,10)
czas=c()
beta=c()
dane=c(5,1,5,14,3,19,1,1,4,22)
czas=c(94.32,15.72,62.88,125.76,5.24,31.44,1.05,1.05,2.10,10.48)
alfa=1.8
beta=c(); beta[1]=0
delta=1
gamma=0.01

for (i in 2:N)
{
	for (j in 1:10)
	lambda[i,j]=rgamma(1,sh=dane[j]+alfa,ra=czas[j]+beta[i-1])
	beta[i]=rgamma(1,sh=gamma+10*alfa,ra=delta+sum(lambda[i,]))
}
par(mfrow=c(1,3))
hist(lambda[,1])
hist(lambda[,2])
hist(beta)

#6.16

#Zastosowanie algorytmu EM do mieszanki rozkładów dla zbioru faithful

library(mixtools)
data(faithful)
attach(faithful)
w=normalmixEM(waiting, arbvar = FALSE, epsilon = 1e-03)
w


#6.17

#Zastosowanie algorytmu EM do mieszanki trzech rozkładów normalnych

library(mixtools)
dane=c(rnorm(20,0,1),rnorm(40,5,1),rnorm(30,12,1))
hist(dane)
w=normalmixEM(dane, arbvar = FALSE, epsilon = 1e-03,k=3)
w