#70.
X<-read.table("C:/Users/Marija/Desktop/OPMS/nerv.txt")
x<-X$V1
n<-length(x)
stat<-function(x){
  return(sqrt(n)*sum((x-mean(x))^3)/(sum((x-mean(x))^2))^(3/2))
}
B<-1000
T<-stat(x)
T #ocena koeficijenta asimetrije
Tb<-c()
for(i in 1:B){
  Xs<-sample(x,n,replace=T)
  Tb[i]<-stat(Xs)
}
t<-mean(Tb)
s2<-sum((Tb-t)^2)/B
se<-sqrt(s2)
se

#intervali poverenja
alfa<-0.05
a<-T-qnorm(alfa/2)*se
b<-T+qnorm(alfa/2)*se
cat("Normalni interval: ", min(a,b),max(a,b))

a1<-sort(Tb)[B*alfa/2]
b1<-sort(Tb)[B*(1-alfa/2)]
cat("Percentil interval: ", a1,b1)

a2<-2*T-sort(Tb)[B*(1-alfa/2)]
b2<-2*T-quantile(Tb)[B*alfa/2]
cat("Stozerni interval: ", a2,b2)

se1<-c()
Tb1<-c()
Tz<-c()
for (b in 1:B) {
  Xs<-sample(x,n,replace=T)
  Tb1[b]<-stat(Xs)
  Xs1<-sample(Xs,B*n,replace = TRUE)
  Xs1<-array(Xs1,c(n,B))
  Tz<-apply(Xs1, 2, stat)
  se1[b]<-sqrt(sum((Tz-mean(Tz))^2)/B)
}
Z<-(Tz-T)/se1
se<-sqrt(sum((Tb-t)^2)/B)
a3<-T-sort(Z)[B*(1-alfa/2)]*se
b3<-T-sort(Z)[B*alfa/2]*se
cat("Studentizovan interval: ", a3,b3)
#drugi nacin za odredjivanje studentizovanog intervala povernja
#preko uticajne krive
sigma2<-1/n*sum((Xs-mean(Xs))^2)
m<-mean(Xs)
Fi<-(Xs-m)^3/(sigma2)^(3/2)-T*(1+3*((Xs-m)^2-sigma2)/(2*sigma2)) #funkcija uticala koef. asimetrije
se1<-sqrt(sum(Fi^2)/n^2)
Z1<-(Tb-T)/se1
a31<-T-quantile(Z1,1-alfa/2)*se
b31<-T-quantile(Z1,alfa/2)*se
cat("Studentizovan interval: ", a31,b31)

#ugradjene funkcije za butstrep
library(boot)
theta.boot <- function(dat,ind) {
  Xs<-dat[ind,1]
  stat(Xs)
}
dim(x)<-c(n,1)
boot.obj<-boot(x, statistic = theta.boot, R = 1000)
print(boot.ci(boot.obj, type = c("basic", "norm", "perc")))