随机模拟2

随机数

随机数检验

  • 安装依赖库
sudo apt-get update && sudo apt-get install -y libgsl-dev libdieharder-dev
  • 安装并加载 RDieHarder 包
if( ! "RDieHarder" %in% list.files(.libPaths())  ) install.packages('RDieHarder')
## Installing package into '/home/travis/R/Library'
## (as 'lib' is unspecified)
library(RDieHarder) # 加载

RDieHarder [@R-RDieHarder] 由 Dirk Eddelbuettel 开发,将 Robert G. Brown 的工作介绍给 R 用户

if( ! "purrr" %in% list.files(.libPaths())  ) install.packages('purrr')
library(purrr)
dieharderGenerators() %>% head
##       names id
## 1  borosh13  0
## 2      cmrg  1
## 3   coveyou  2
## 4 fishman18  3
## 5 fishman20  4
## 6 fishman2x  5
dieharderTests() %>% head
##                names id
## 1  diehard_birthdays  0
## 2     diehard_operm5  1
## 3 diehard_rank_32x32  2
## 4   diehard_rank_6x8  3
## 5  diehard_bitstream  4
## 6       diehard_opso  5
set.seed(2018) 
n <- 2^24
x <- runif(n,0,1)
delta <- 0.01
len <- diff(c(0,which(x < delta),n+1))-1
ylim <- seq( 0, 1800, by = 300)
xlim <- seq( 0, 100, by = 20)
p <- hist(len[len < 101], breaks = -1:100+0.5, plot = FALSE)
par(mar = c(2,2,.5,.5))
plot(p, xlab = '间距', ylab = '频数', axes = FALSE, 
     col = "lightblue", border = "white", main = "")     
axis( 1, labels = xlim, at = xlim, las = 1) # x 轴
axis( 2, labels = ylim, at = ylim, las = 0) # y 轴
box(col="gray")
# 添加线性回归线
xx <- seq.int(from = 0, to = 100, by = 1)
xy <- p$counts
options(digits = 2)
fit <- lm(xy~xx)
abline(fit, col = 'red')
游程直方图^[选自:[随机数生成及其在统计模拟中的应用](https://cosx.org/2017/05/random-number-generation/)]

Figure 1: 游程直方图1

b = coef(fit)
# mtext(paste0( "Y = ", paste0(paste0(b[1], b[2]),"x") ), side = 3, cex = 2)

并行模拟

如何解决并行模拟的问题

  • parallel

  • Rmpi

  • MCMC 并行

  • Stan

假设有4个人,每个人手里有一张100块,每次都随机地给他人1块,若手里没钱了,就不给了,最后钱的分布

\[P = \begin{bmatrix}0 & \frac{1}{3} & \frac{1}{3} & \frac{1}{3}\\ \frac{1}{3} & 0 & \frac{1}{3} & \frac{1}{3}\\ \frac{1}{3} & \frac{1}{3} & 0 & \frac{1}{3} \end{bmatrix}\]

# n 人数
n <- 10
# n 阶单位矩阵
I_mat <- diag(rep(1, n))
# n 阶全1矩阵
One_mat <- rep(1, n) %o% rep(1, n)
# 概率转移矩阵
P_mat <- 1 / (n - 1) * (One_mat - I_mat) 

布朗运动

# code from https://yihui.name/cn/2009/02/cumsum-of-normal-var-and-uniform-var/
n = 1e+05
x = cumsum(rnorm(n, 0, 1))
y = cumsum(runif(n, -1, 1))
cl = gray(seq(0, 1, length = n))
plot(x, y, pch = ".", col = cl, xlab = 'CUM-Norm',
  ylab = 'CUM-Uniform')

样条曲线

x <- sort(stats::rnorm(5))
y <- sort(stats::rnorm(5))
plot(x, y, pch = 19)
res <- xspline(x, y, 1, draw = TRUE) #  X-spline
lines(res)

相关

comments powered by Disqus