###for calculating running mean
runningmean = function (x,N){
  #y = rpois(N,x)
  y =runif(N,2,3)
  
  c = cumsum(y)
  n = 1:N
  c/n
}
##plotting running mean and Simulating Law of Large Numbers
library(viridis)
par(mfrow=c(1,3))
u = runningmean(1, 100)
x=1:100; plot(u~x, type="l");
replicate(10, 
          lines(runningmean(1, 100)~x, type="l", col = sample(viridis(15, option="A",1))
          )
)
u = runningmean(1, 1000)
x=1:1000; plot(u~x, type="l");
replicate(10, lines(runningmean(1, 1000)~x, type="l",  col = sample(viridis(15, option="B",1))
)
)
u = runningmean(1, 10000)
x=1:10000; plot(u~x, type="l");
replicate(10, lines(runningmean(1, 10000)~x, type="l",   col = sample(viridis(15, option="C",1))
)
)

##Simulating Binomial Central Limit Theorem
binomialsim1 = rbinom(100,10,0.1)
# generates 100 Binomial (10,0.1) samples
binomialsim2 = rbinom(100,10,0.25)
# generates 100 Binomial (10,0.25) samples
binomialsim3 = rbinom(100,10,0.5)
# generates 100 Binomial (10,0.5) samples

hist(binomialsim1, main = "Binomial(10, 0.1)")
hist(binomialsim2, main = "Binomial(10, 0.25)")
hist(binomialsim3, main = "Binomial(10, 0.5)")


binomialsim1 = rbinom(100,10,0.1)
# generates 100 Binomial (10,0.1) samples
stbinomialsim1= (binomialsim1-10*0.1)/(sqrt(10*0.1*0.9))

binomialsim1a = rbinom(100,100,0.1)
# generates 100 Binomial (100,0.1) samples
stbinomialsim1a= (binomialsim1a-100*0.1)/(sqrt(100*0.1*0.9))


binomialsim1b = rbinom(100,1000,0.1)
# generates 100 Binomial (1000,0.1) samples
stbinomialsim1b= (binomialsim1b-1000*0.1)/(sqrt(1000*0.1*0.9))


par(mfrow=c(1,3))
hist(stbinomialsim1, main = "std-Binomial(10, 0.1)")
hist(stbinomialsim1a, main = "std-Binomial(100, 0.1)")
hist(stbinomialsim1b, main = "std-Binomial(1000, 0.1)")

x = rnorm(100)
par(mfrow=c(1,2))
qqnorm(x)
qqline(x)
qqnorm(stbinomialsim1b)
qqline(stbinomialsim1b)

par(mfrow=c(1,3))
boxplot(x, stbinomialsim1)
boxplot(x, stbinomialsim1a)
boxplot(x, stbinomialsim1b)

##towards Simulating Central Limit Theorem
STDS1000= (binomialsim1b-1000*0.1)/(sqrt(1000*0.1*0.9))
par(mfrow=c(1,3))
qqnorm(STDS1000)
qqline(STDS1000)
boxplot(x, STDS1000)
hist(STDS1000, main = "STD-1000")

### CLT For uniform
Mean = replicate(100, mean(runif(100, 10)))
STDMean = sqrt(100)*(Mean -0.1)/(sqrt(0.01))
hist(STDMean, main = "Standardised Uniform", xlab="x")
x = rnorm(100)
par(mfrow=c(1,3))
boxplot(x,STDMean)
qqnorm(x);qqline(x)
qqnorm(STDMean);qqline(STDMean)


### CLT For exponential
Mean = replicate(100, mean(runif(100, 10)))
STDMean = sqrt(100)*(Mean -0.1)/(sqrt(0.01))
hist(STDMean, main = "Standardised Uniform", xlab="x")
x = rnorm(100)
par(mfrow=c(1,3))
boxplot(x,STDMean)
qqnorm(x);qqline(x)
qqnorm(STDMean);qqline(STDMean)
library(moments)
c(skewness(x),skewness(STDMean))
c(kurtosis(x), kurtosis(STDMean))