#
#   https://thalis.math.upatras.gr/~vpiperig/
#           Math-Biost/Distributions_Additional_2.pdf
#
#   Kwdikas gia tnv prosomoiwsn kai ta sxnmata: sel. 14, 16
#
#

#tuxaio deigma 1000 paratnrnsewv N(0,1)
x1<-rnorm(1000)

# istogramma tns X1^2 (embadov=1) 
hist(x1^2,nclass=25,col=4,prob=TRUE)

# pros8nkn tns suvartnsns pukvotntas tns x ^2 katavomns
curve(dchisq(x,df=1),lwd=3,col="red",add=TRUE)

x2<-rnorm(1000)
hist(x1^2+x2^2,nclass=25,col=4,prob=TRUE)
curve(dchisq(x,df=2),lwd=3,col="red",add=TRUE)

x3<-rnorm(1000)
x4<-rnorm(1000)
x5<-rnorm(1000)
hist(x1^2+x2^2+x3^2+x4^2+x5^2,nclass=25,col=4,prob=TRUE)
curve(dchisq(x,df=5),lwd=3,col="red",add=TRUE)

# oi deigmatikes meses times kai oi diaspores, briskontai polu kovta 
# stis 8ewpitikes times E(X^2_n)=n kai Var(X^2_n)=2n
# edw n=1 n=2 n=3 n=5  
round(c(mean(x1^2),mean(x1^2+x2^2),mean(x1^2+x2^2+x3^2),mean(x1^2+x2^2+x3^2+x4^2+x5^2)),2)
round(c(var(x1^2),var(x1^2+x2^2),var(x1^2+x2^2+x3^2),var(x1^2+x2^2+x3^2+x4^2+x5^2)),2)

################################################################
#
#  pairvoume 5 deigmata mege8ous n=3

X1<-rnorm(5)
X2<-rnorm(5)
X3<-rnorm(5)

rbind(X1,X2,X3)

# upologizetai n deigmatikn mesn timn gia ka8va apo ta 5 deigmata 
m<-apply(rbind(X1,X2,X3),2,mean)
m
# upologizetai n deigmatikn diaspora gia ka8va apo ta 5 deigmata 
S2<-apply(rbind(X1,X2,X3),2,var)
S2
##################################################################
#
#  pairvoume 1000 deigmata (gia tnv kalutern patnrnsn tns 8ewpias)
#  mege8ous n=4 apo N(0,1)

X1<-rnorm(1000)
X2<-rnorm(1000)
X3<-rnorm(1000)
X4<-rnorm(1000)

m<-apply(rbind(X1,X2,X3,X4),2,mean)
S2<-apply(rbind(X1,X2,X3,X4),2,var)

# (n-1)S^2/sigma^2=3*S^2/1
#
hist((4-1)*S2/1^2,nclass=50,col=4,prob=TRUE)
curve(dchisq(x,df=3),lwd=3,col="red",add=TRUE)
# (bar(X)-mu)/(S/sqrt(n))=(bar(X)-0)/(S/sqrt(4))
#
hist((m-0)/sqrt(S2/4),nclass=50,col=4,prob=TRUE)
curve(dt(x,df=3),lwd=3,col="red",add=TRUE)

#######################################################################
#  pairvoume 1000 deigmata 
#  mege8ous n=4 apo N(-1,2^2)

X1<-rnorm(1000,mean=-1,sd=2)
X2<-rnorm(1000,mean=-1,sd=2)
X3<-rnorm(1000,mean=-1,sd=2)
X4<-rnorm(1000,mean=-1,sd=2)

m<-apply(rbind(X1,X2,X3,X4),2,mean)
S2<-apply(rbind(X1,X2,X3,X4),2,var)

# (n-1)S^2/sigma^2=3*S^2/4
#
hist((4-1)*S2/2^2,nclass=50,col=4,prob=TRUE)
curve(dchisq(x,df=3),lwd=3,col="red",add=TRUE)

# (bar(X)-mu)/(S/sqrt(n))=(bar(X)-(-1))/(S/sqrt(4))
hist((m-(-1))/sqrt(S2/4),nclass=50,col=4,prob=TRUE)
curve(dt(x,df=3),lwd=3,col="red",add=TRUE)

#######################################################################