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)

integrate(f2, 0, Inf)
## 1 with absolute error < 5.5e-05

Площадь под кривой равна 1

h1 <- hist(v, freq = FALSE)

lines(x2,f2(x2), col="blue", lwd=2)

MCMC

Norm

s <- 0
x=rep(runif(1),Nsim) # initialize the chain
for (i in 2:Nsim){
      
      ## norm
      y = rnorm(1)
      rho = f2(y)*dnorm(x[i-1])/(f2(x[i-1])*dnorm(y))
      
      tmp <- runif(1)
      x[i]=x[i-1] + (y-x[i-1])*(tmp<rho)
      s <- s + (tmp<rho)
}

burnIn <- 2000
h2 <- 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
lines(x2, f2(x2), col="blue", lwd=2) 

acceptance_rate <- s/Nsim
print(acceptance_rate)
## [1] 0.09385714
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.7784, p-value < 2.2e-16
## alternative hypothesis: two-sided

Cauchy

s <- 0
x=rep(runif(1),Nsim) # initialize the chain
for (i in 2:Nsim){
      
      ## cacuhy
      y = rt(1, 1)
      rho = f2(y)*dt(x[i-1], 1)/(f2(x[i-1])*dt(y, 1))
      
      tmp <- runif(1)
      x[i]=x[i-1] + (y-x[i-1])*(tmp<rho)
      s <- s + (tmp<rho)
}

burnIn <- 5000
h2 <- 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
lines(x2, f2(x2), col="blue", lwd=2) 

acceptance_rate <- s/Nsim
print(acceptance_rate)
## [1] 0.1491429
ks.test(jitter(x[-c(1:burnIn)]),jitter(v))
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  jitter(x[-c(1:burnIn)]) and jitter(v)
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided

Exp

s <- 0
x=rep(runif(1),Nsim)
for (i in 2:Nsim){
      
      ## exp
      y = rexp(1)
      rho = f2(y)*dexp(x[i-1])/(f2(x[i-1])*dexp(y))
      
      tmp <- runif(1)
      x[i]=x[i-1] + (y-x[i-1])*(tmp<rho)
      s <- s + (tmp<rho)
}

burnIn <- 5000
h2 <- 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
lines(x2, f2(x2), col="blue", lwd=2) 

acceptance_rate <- s/Nsim
print(acceptance_rate)
## [1] 0.3307143
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.04, p-value = 0.01378
## alternative hypothesis: two-sided

Gamma

s <- 0
x=rep(runif(1),Nsim)
for (i in 2:Nsim){
      
      ## gamma
      y = rgamma(1,1)
      rho = f2(y)*dgamma(x[i-1],1)/(f2(x[i-1])*dgamma(y,1))
      
      tmp <- runif(1)
      x[i]=x[i-1] + (y-x[i-1])*(tmp<rho)
      s <- s + (tmp<rho)
}

burnIn <- 5000
h2 <- 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
lines(x2, f2(x2), col="blue", lwd=2) 

acceptance_rate <- s/Nsim
print(acceptance_rate)
## [1] 0.337
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.058214, p-value = 5.274e-05
## alternative hypothesis: two-sided