image-20230408163433029

算法：

生成一个[0,1][0,1]上的均匀分布的随机数，再将[0,1][0,1]划分为两个区间，判断其落入哪个区间，确定其取值。

根据中心极限定理，当nn足够大时，Yn∼N(np,np(1−p))Y_n \sim N(np, np(1-p))，其中p=1/3p=1/3，Yn/n=fnY_n/n=f_n，对fnf_n进行标准化，得

Zn=(fn−p)p(1−p)/n=(Yn−np)np(1−p)Z_n=(f_n-p)\sqrt {p(1-p)/n}=(Y_n-np)\sqrt{np(1-p)}

ZnZ_n服从标准正态分布，故当nn充分大时，fnf_n的渐进分布为N∼(13,29n)N \sim (\frac{1}{3},\sqrt{\frac{2}{9n}})

R程序：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
rm(list=ls())
rng <- function(n){
  xvec <- numeric(n)  # 储存结果的向量
  for(i in 1:n){
    U <- runif(1)
    if(U <= 1/3){  # 判断U落入区域
      xvec[i] <- 1
    }else{
      xvec[i] <- 2
    }
  }
  xvec
}
# 生成随机数
x1 <- rng(100)
x2 <- rng(1000)
x3 <- rng(10000)
# x取1的百分比
N1 <- sum(x1 == 1)/100
N2 <- sum(x2 == 1)/1000
N3 <- sum(x3 == 1)/10000
paste("n = 100, N = ", N1)
paste("n = 1000, N = ", N2)
paste("n = 10000, N = ", N3)

## 检验是否服从渐进分布
# 定义函数，计算fn
fn <- function(n, x){
  result <- rng(n)
  sum(result == x)/n
}
# 定义函数，检验是否服从渐进分布
Central_Limit_Theorem <- function(N, n){
  f_n <- replicate(N, fn(n, 1))  # 调用函数N次
  ks.test(f_n, "pnorm", mean = 1/3, sd = sqrt(2/(9*n)))
}
# 
Central_Limit_Theorem(100, 100)
Central_Limit_Theorem(100, 1000)
Central_Limit_Theorem(100, 10000)

结果：

n=100,1000,10000n=100,1000,10000的随机数X取1的百分比分别为

image-20230408183852627

对于n=100,1000,10000n=100,1000,10000，分别生成100组数据，并计算其百分比fnf_n，再利用K-S检验，检验其是否符合渐进分布，其检验结果如下

（1）n=100n=100

image-20230408184213162

（2）n=1000n=1000

image-20230408184236169

（3）n=10000n=10000

image-20230408184309857

image-20230407165258305

1.对分布函数p(x)=1nx1n−1,x∈[0,1]\begin{aligned} p(x) = \frac{1}{n} x^{\frac{1}{n} - 1}, x \in [0,1] \end{aligned}求积分，得F(x)=x1nF(x)=x^{\frac{1}{n}}，故反函数F−1(y)=ynF^{-1}(y)=y^{n}

R程序：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 第一问
rm(list=ls())
set.seed(0721) # 设置随机数种子
# 定义参数n
n <- 3
# 定义样本大小
N <- 10000
rng.F1 <- function(N, n){
  runif(N)^n  # 反函数
}
x1 <- rng.F1(N, n)
# 绘制直方图和理论密度曲线
hist(x1, prob = TRUE, main = expression(paste("Beta(", frac(1,n), ", 1)", " distribution")))
curve(dbeta(x, 1/n, 1), add = TRUE, col = "red")

结果：

取n=3n=3，模拟生成10000个服从Beta(1n,1)Beta(\frac{1}{n},1)分布的随机数，直方图如下

image-20230407162743519

2.对分布函数p(x)=nxn−1,x∈[0,1]\begin{aligned} p(x) = n x^{n - 1}, x \in [0,1] \end{aligned}，积分得F(x)=xnF(x)=x^{n}，故反函数F−1(y)=y1nF^{-1}(y)=y^{\frac{1}{n}}

R程序：

1
2
3
4
5
6
7
8
9
10
11
12
# 第二问
# 定义参数n
n <- 3
# 定义样本大小
N <- 10000
rng.F2 <- function(N, n){
  runif(N)^(1/n)  # 反函数
}
x2 <- rng.F2(N, n)
# 绘制直方图和理论密度曲线
hist(x2, prob = TRUE, main = expression(paste("Beta(n, 1)", " distribution")))
curve(dbeta(x, n, 1), add = TRUE, col = "red")

结果：

取n=3n=3，模拟生成10000个服从Beta(n,1)Beta(n,1)分布的随机数，直方图如下

image-20230407163025997

3.对分布函数p(x)=2π1−x2,x∈[0,1]\begin{aligned} p(x) = \frac{2}{\pi \sqrt{1 - x^2}}, x \in [0,1] \end{aligned}，积分得F(x)=2πarcsinxF(x)=\frac{2}{\pi}arcsinx，故反函数F−1(y)=sin(π2y)F^{-1}(y)=sin(\frac{\pi}{2}y)

R程序：

1
2
3
4
5
6
7
8
9
10
# 第三问
# 定义样本大小
N <- 10000
rng.F3 <- function(N){
  sin( pi/2*runif(N) )  # 反函数
}
x3 <- rng.F3(N)
# 绘制直方图和理论密度曲线
hist(x3, prob = TRUE, main = expression(paste(p(x))))
curve(2/(pi*sqrt(1-x^2)), add = TRUE, col = "red")

结果：

模拟生成10000个服从p(x)p(x)的随机数，直方图如下

image-20230407163447793

4.对分布函数p(x)=1π(1+x2),x∈(−∞,∞)\begin{aligned} p(x) = \frac{1}{\pi (1 + x^2)}, x \in (-\infty, \infty) \end{aligned}，积分得F(x)=1πarctanxF(x)=\frac{1}{\pi}arctanx，故反函数F−1(y)=sin(πy)F^{-1}(y)=sin(\pi y)

R程序：

1
2
3
4
5
6
7
8
9
# 第四问
# 定义样本大小
N <- 10000
rng.F4 <- function(N){
  tan( pi*runif(N) )  # 反函数
}
x4 <- rng.F4(N)
hist(x4, prob = TRUE, main = expression(paste("Cauchy distribution")))
curve(1/(pi*(1+x^2)), add = TRUE, col = "red")

结果：

模拟生成10000个服从柯西分布的随机数，直方图如下

image-20230407163848566

5.对分布函数p(x)=cos⁡(x),x∈[0,π2]\begin{aligned} p(x) = \cos(x), x \in [0, \frac{\pi}{2}] \end{aligned}，积分得F(x)=sinxF(x)=sinx，故反函数F−1(y)=arcsinyF^{-1}(y)=arcsiny

R程序：

1
2
3
4
5
6
7
8
9
10
# 第五问
# 定义样本大小
N <- 10000
rng.F5 <- function(N){
  asin( runif(N) )  # 反函数
}
x5 <- rng.F5(N)
# 绘制直方图和理论密度曲线
hist(x5, prob = TRUE, main = expression(paste(p(x))))
curve(cos(x), add = TRUE, col = "red")

结果：

模拟生成10000个服从p(x)p(x)的随机数，直方图如下

image-20230407164110699

6.对分布函数p(x)=αηxα−1e−xαη, x>0 (α>0,η>0)\begin{aligned} p(x) = \frac{\alpha}{\eta} x^{\alpha-1} e^{-\frac{x^\alpha}{\eta}}, \ x>0 \ (\alpha>0, \eta>0) \end{aligned}，积分得F(x)=1−e−xαηF(x)=1-e^{-\frac{x^\alpha}{\eta}}，故反函数F−1(y)=[−ηln(1−y)]1αF^{-1}(y)=[-\eta ln(1-y)]^{\frac{1}{\alpha}}

R程序：

1
2
3
4
5
6
7
8
9
10
11
12
# 第六问
# 定义参数
alpha <- 2
lambda <- 1
# 定义样本大小
N <- 10000
rng.F6 <- function(N, alpha, lambda){
  ( (-1)*lambda*log(1 - runif(N)) )^(1/alpha)  # 反函数
}
x6 <- rng.F6(N, alpha, lambda)
hist(x6, prob = TRUE, main = expression(paste("Weibull distribution")))
curve((alpha/lambda)*(x^(alpha-1))*exp((-1)*x^(alpha)/lambda), add = TRUE, col = "red")

结果：

取α=2,η=1\alpha =2,\eta =1，模拟生成10000个服从分威布尔布的随机数，直方图如下

image-20230407165059945

image-20230407165352675

算法：

将[0,1][0,1]区间划分为10段，生成一个服从[0,1][0,1]上的均匀分布的随机数，判断其落入哪个区间，确定其取值

R程序：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
rm(list=ls())
set.seed(0721) # 设置随机数种子, Ciallo~
# 输入为目标分布的取值和概率
rng <- function(x, prob) {
  # 计算目标分布的累积分布函数
  F <- cumsum(prob)
  U1 <- runif(1)
  # 根据U1的值确定目标分布属于哪一个离散均匀分布函数
  i <- which(F > U1)[1]
  X <- x[i]
  return(X)
}

# 生成10000个随机数
x <- 1:10 # 目标分布的取值
prob <- c(0.06, 0.06, 0.06, 0.06, 0.06, 0.15, 0.13, 0.14, 0.15, 0.13) # 目标分布的概率
result <- replicate(10000, rng(x, prob))  # 调用函数10000次
table(result)

结果：

image-20230407170706974

模拟生成10000个随机数，可以看出，各个取值频率基本与概率一致

image-20230408150155141

算法：

1.逆变换法

对分布密度函数p(x)=12e−∣x∣, −∞<x<∞\begin{aligned} p(x) = \frac12 e^{-|x|}, \ -\infty < x < \infty \end{aligned}积分，得到累积密度函数

F(x)={1−12e−x,x≥012ex,x<0\begin{aligned} F(x) =& \begin{cases} 1-\frac{1}{2}e^{-x}, & x \geq 0 \\ \frac{1}{2}e^{x}, & x < 0 \end{cases} \end{aligned}

则反函数为

F−1(u)={−ln[2(1−u)],u≥0ln(2u),u<0.\begin{aligned} F^{-1}(u) = \begin{cases} -ln[2(1-u)], & u \geq 0 \\ ln(2u), & u < 0 . \end{cases} \end{aligned}

2.复合法

将p(x)p(x)拆分为g1(x)=e−x,x≥0g_1(x)=e^{-x},x \geq 0，g2(x)=ex,x<0g_2(x)=e^x,x < 0，即p(x)=0.5∗g1(x1)+0.5∗g2(x2)p(x)=0.5*g_1(x_1)+0.5*g_2(x_2)。又通过逆变换法，求得g1−1(x)=−ln(2y)g^{-1}_1(x)=-ln(2y)，g2−1(x)=ln(2y)g^{-1}_2(x)=ln(2y)，通此函数可生成g1(x)g_1(x)，g2(x)g_2(x)的随机数，再设置随机数，依概率生成即可

R程序：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
rm(list=ls())
set.seed(0721) # 设置随机数种子
# 逆变换法
N <- 10000
rng.F1 <- function(N){
  U <- runif(N)
  ifelse(U >= 0.5, (-1)*log(2*(1 - U)), log(2*U))
}
x1 <- rng.F1(N)
# 绘制直方图和理论密度曲线
hist(x1, prob = TRUE, main = expression(paste("Laplace distribution inverse transform method")))
curve(exp(-abs(x))/2, add = TRUE, col = "red")

# 复合法
rng.F2 <- function(n){
  xvec <- numeric(n)  # 储存结果的向量
  for(i in 1:n){
    
    U <- runif(1)
    if(U <= 0.5){  # 判断U落入区域
      xvec[i] <- -log(2*runif(1))
    }else{
      xvec[i] <- log(2*runif(1))
    }
  }
  xvec
}
x2 <- rng.F2(N)
# 绘制直方图和理论密度曲线
hist(x2, prob = TRUE, main = expression(paste("Laplace distribution Composition method")))
curve(exp(-abs(x))/2, add = TRUE, col = "red")

结果：

1.逆变换法

image-20230408155257808

2.复合法

image-20230408155325203

image-20230408155622313

算法：

采用舍选法，根据分析，p(x)p(x)在定义域内单调递减，故最大值M=p(0,8)M=p(0,8)，再根据舍选法I的算法，如下图

image-20230408163319434

即可模拟得随机数

R程序：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
rm(list=ls())
set.seed(0721) # 设置随机数种子
# 舍选法
N <- 10000
rng.F1 <- function(n){
  xvec <- numeric(n)  # 存储结果的向量
  for(i in 1:n){
    repeat{
      U1 <- runif(1)
      U2 <- runif(1)
      X <- 0.8 + 0.2*U1
      Y <- (1/0.000336)*0.8*0.2^3*U2  # M为3/2
      if(Y <= 1/0.000336*X*(1 - X)^3) break  # 判断是否接受
    }
    xvec[i] <- X
  }
  
  xvec
}
x1 <- rng.F1(N)
# 绘制直方图和理论密度曲线
hist(x1, prob = TRUE, main = expression(paste(p(x))))
curve(1/0.000336*x*(1 - x)^3, add = TRUE, col = "red")

结果：

image-20230408163353606