#1. butstrep za medijanu
X<-rnorm(10)
T<-median(X)
Tb<-c()
for(i in 1:50){
  Xs<-sample(X,50,replace=T)
  Tb[i]<-median(Xs)
}
t<-mean(Tb)
s2<-sum((Tb-t)^2)/50
se<-sqrt(s2)
se

#2. Izrachunati 10 vrednosti ocene lambda_n(Fn) za normalne mesavine
#oblika F(x)=(1-epsilon)Fi(x)+epsilon*Fi(x/tau), za 4 kombinacije 
#epsilon=0.1,0.2, tau=2,3, uzimanjem uzorka n=50 iz date raspodele i
#onda uzimanjem 10 butstrep uzorka velicine B=100 iz Fn.

n<-50
B<-100
#uzorak iz normalne mesavine
X<-c()
komponente<-sample(1:2,prob=c(0.1,0.9),size=1000,replace=TRUE)
m<-c(0,0)
sd<-c(1,2)
X<-rnorm(n,m[komponente],sd[komponente]) #uzorak iz normalne mesavine 0.1*Fi(x)+0.9*Fi(x/2)

T<-matrix(0,10,B) 
Tb<-c()
Xs<-matrix(0,B,n) #matrica u koju se smestaju uzorci
for(j in 1:10){
  for(i in 1:B){
    Xs[i,]<-sample(X,n,replace=TRUE)
    Tb[i]<-(sqrt(n)*(mean(Xs[i,])-mean(X)))/sqrt((1/n*sum(X-mean(X))^2)) 
    T[j,i]<-Tb[i]
  }
}
#Xs se moze napraviti kao matrica velicine n*B, a zatim se mogu 
#koristiti funkcije rowSum ili rowMeans za racunanje potrebnih 
#vrednosti kako se ne bi koristile dve for petlje

a<-0 #zadata vrednost
lambda<-rep(0,10) #matrica u koju se smesta 10 ocena 
for(j in 1:10){
  for(i in 1:B){
    if(T[j,i]<=a) #ako je zadovoljen uslov u okviru verovatnoce, povecava se frekvencija
      lambda[j]=lambda[j]+1/B
    
  }
}
lambda

#3.Reshiti prethodni problem kada 10 vrednosti lambda_n(Fn) su 
#dobijene uzimanjem novog uzorka (X_i1,...,X_in), i=1,...,10, za 
#svaki od 10 slucajeva i jednog butstrepa B=100 za svaki uzorak
#x_i1,...,x_in.

n<-50
B<-100
#uzorak iz normalne mesavine 0.1*Fi(x)+0.9*Fi(x/2)
X<-matrix(0,10,n)
komponente<-sample(1:2,prob=c(0.1,0.9),size=1000,replace=TRUE)
m<-c(0,0)
sd<-c(1,2)
for(i in 1:10){
  X[i,]<-rnorm(n,m[komponente],sd[komponente]) 
}

T<-matrix(0,10,B) 
Tb<-c()
Xs<-matrix(0,B,n) #matrica u koju se smestaju uzorci
for(j in 1:10){
  for(i in 1:B){
    Xs[i,]<-sample(X[j,],n,replace=TRUE)
    Tb[i]<-(sqrt(n)*(mean(Xs[i,])-mean(X[j,])))/sqrt((1/n*sum(X[j,]-mean(X[j,]))^2)) 
    T[j,i]<-Tb[i]
  }
}

a<-0 #zadata vrednost
lambda<-rep(0,10) #matrica u koju se smesta 10 ocena 
for(j in 1:10){
  for(i in 1:B){
    if(T[j,i]<=a) #ako je zadovoljen uslov u okviru verovatnoce, povecava se frekvencija
      lambda[j]=lambda[j]+1/B
    
  }
}
lambda

#4. Odrediti ocenu koeficijenta asimetrije, ocenu standardne 
#devijacije ocene i intervale poverenja za podatke iz baze nerv

X<-read.table("C:/Users/Marija/Desktop/OPS/nerv.txt")
x<-c(X$V1,X$V2,X$V3,X$V4,X$V5,X$V6)
x1<-x[!is.na(x)]
n<-length(x1)
#funkcija za rachunanje koeficijenta asimetrije
stat<-function(x){
  return(sqrt(n)*sum((x-mean(x))^3)/(sum((x-mean(x))^2))^(3/2))
}
T<-stat(x1)
T #ocena koeficijenta asimetrije
B<-1000
Tb<-c()
for(i in 1:B){
  Xs<-sample(x1,n,replace=T)
  Tb[i]<-stat(Xs)
}
t<-mean(Tb)
s2<-sum((Tb-t)^2)/B
se<-sqrt(s2)
se #ocena standardne devijacije

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

a1<-quantile(Tb,alfa/2)
b1<-quantile(Tb,1-alfa/2)
cat("Percentil interval: ", a1,b1)

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

se1<-c()
Tb<-c()
Tz<-c()

for (i in 1:1000) {
  Xs<-sample(x1,n,replace=T)
  Tb[i]<-stat(Xs)
  for(j in 1:1000){
    Xs1<-sample(Xs,n,replace=T)
    Tz[j]<-stat(Xs1)
  }
  se2<-sqrt(sum((Tz-mean(Tz))^2)/B)
  se1[i]<-se2
}
Z<-(Tb-T)/se1
se<-sqrt(sum((Tb-t)^2)/B)
a3<-T-quantile(Z,1-alfa/2)*se
b3<-T-quantile(Z,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(x1)<-c(n,1)
boot.obj<-boot(x1, statistic = theta.boot, R = 1000)
print(boot.ci(boot.obj, type = c("basic", "norm", "perc")))