MCMC_kolmogorov_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)

metropolis_hasting <- function(proposal, Nsim, param){
      
      #proposal function
      if (proposal == "n") {proposalfun =  function(x, param = c(0,1)) return(dnorm(x, param[1], param[2]))
      rproposalfun =  function(param = c(0,1)) return(rnorm(1, param[1], param[2]))}
      
      if (proposal == "u") {proposalfun =  function(x, param = c(0,1)) return(dunif(x, param[1], param[2]))
      rproposalfun =  function(param = c(0,1)) return(runif(1, param[1], param[2]))}
      
      if (proposal == "b") {proposalfun =  function(x, param = c(1,1)) return(dbeta(x,param[1],param[2]))
      rproposalfun =  function(param = c(1,1)) return(rbeta(1, param[1], param[2]))}
      
      if (proposal == "g") {proposalfun =  function(x, param = c(1,1)) return(dgamma(x,param[1], param[2]))
      rproposalfun =  function(param = c(1,1)) return(rgamma(1, param[1], param[2]))}
      
      if (proposal == "c") {proposalfun =  function(x, param = c(0,1)) return(dcauchy(x, param[1], param[2]))
      rproposalfun =  function(param = c(0,1)) return(rcauchy(1, param[1], param[2]))}
      
      if (proposal == "e") {proposalfun =  function(x, param = 1) return(dexp(x,param[1]))
      rproposalfun =  function(param = 1) return(rexp(1, param[1]))}
      
      #initial values
      s1 <- 0
      x=rep(runif(1),Nsim) # initialize the chain
      
      f2x <- f2(x[1])
      y=rproposalfun(param)
      f2y <- f2(y)
      
      rho=f2y*proposalfun(x[1], param)/(f2x*proposalfun(y, param))
      tmp = runif(1)
      
      
      for (i in 2:Nsim){
            if (tmp<rho) {
                  f2x <- f2y
                  y=rproposalfun(param)
                  f2y <- f2(y)
                  rho=f2y*proposalfun(x[i-1], param)/(f2x*proposalfun(y, param)) 
                  
                  tmp <- runif(1)
                  x[i]=x[i-1] + (y-x[i-1])*(tmp<rho)
                  s1 <- s1 + 1
            }
            else
            {
                  y=rproposalfun(param)
                  f2y <- f2(y)
                  rho=f2y*proposalfun(x[i-1], param)/(f2x*proposalfun(y, param)) 
                  
                  tmp <- runif(1)
                  x[i]=x[i-1] + (y-x[i-1])*(tmp<rho)
            }
      }
      return(x)
}

x <- metropolis_hasting("c", 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, f2(x2), 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.030943, p-value = 0.007506
## 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.026771, p-value = 0.03057
## alternative hypothesis: two-sided