MCMC_Kol_norm

Density

Выборка из распределения Колмогорова, посчитанная напрямую

Nsim <- 7000

n <- 4000
v <- rep(0,n)
for (i in 1:Nsim){
      x1 <- rnorm(n)
      v[i] <- ks.test(x1, "pnorm")$statistic*sqrt(n)
}

Формула плотности

f2 <- function(x, n = 7000){
      res = 0
      for (k in 1:n){
            res = res + 8*(-1)^(k-1)*k^2*x*exp(-2*k^2*x^2)
      }
      return(res)
}

x2 <- seq(0,2, by = 0.01)

h1 <- hist(v, freq = FALSE)
y2 <- f2(x2)
lines(x2,y2, col="blue", lwd=2)

x <- metropolis_hasting("n", Nsim, param = c(0,1))
x_ <- metropolis_hasting("c", Nsim, param = c(0.85,1))

burnIn <- 2000
h2 <- hist(x[-c(1:burnIn)], plot = FALSE)
h3 <- hist(x_[-c(1:burnIn)], plot = FALSE)
plot( h1, freq = FALSE, col=rgb(0,0,1,1/4), xlim=c(0,2), ylim = c(0,2.5))  # first histogram
plot( h2, freq = FALSE, col=rgb(1,0,0,1/4), xlim=c(0,2), add = T)  # second
plot(h3, freq = FALSE, col=rgb(0,1,0,1/4), xlim=c(0,2), add = T)
lines(x2, y2, col="blue", lwd=2) 

#acceptance_rate <- s/Nsim
#print(acceptance_rate)

ks.test(jitter(x[-c(1:burnIn)]),jitter(v))

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  jitter(x[-c(1:burnIn)]) and jitter(v)
## D = 0.0482, p-value = 2.602e-06
## alternative hypothesis: two-sided

ks.test(jitter(x_[-c(1:burnIn)]),jitter(v))

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  jitter(x_[-c(1:burnIn)]) and jitter(v)
## D = 0.5352, p-value < 2.2e-16
## alternative hypothesis: two-sided