Density

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

Nsim <- 10000

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 = 4000){
      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

Unif

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


burnIn <- 5000
h1 <- hist(v, freq = FALSE)

h2 <- hist(x[-c(1:burnIn)], freq = 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.4766
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.2676, p-value < 2.2e-16
## alternative hypothesis: two-sided

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 <- 5000
h1 <- hist(v, freq = FALSE)

h2 <- hist(x[-c(1:burnIn)], freq = 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.069
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.4408, 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
h1 <- hist(v, freq = FALSE)

h2 <- hist(x[-c(1:burnIn)], freq = 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.142
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.0301, p-value = 0.004763
## 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
h1 <- hist(v, freq = FALSE)

h2 <- hist(x[-c(1:burnIn)], freq = 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.3364
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.0176, p-value = 0.2531
## alternative hypothesis: two-sided

Beta

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

burnIn <- 5000
h1 <- hist(v, freq = FALSE)

h2 <- hist(x[-c(1:burnIn)], freq = 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.4663
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.2678, p-value < 2.2e-16
## 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
h1 <- hist(v, freq = FALSE)

h2 <- hist(x[-c(1:burnIn)], freq = 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.3325
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.0176, p-value = 0.2531
## alternative hypothesis: two-sided