R语言：作业六（逆变换法生成随机变量；线性同余发生器LCG的编写）

news/2024/10/19 4:21:59/

1.

用逆变换法编写产生下述随机变量的程序：
X 0 1
p 0.4 0.6
模拟10000次，并确定随机变量的值0的比例。

f <- function(n){{x <- rep(0,n)for (i in 1:n){u <- runif(1)if(u < 0.6) x[i] <- 1else x[i] <- 0}x}
}
test <- f(10000)
table(test)

test
0 1
4087 5913

a <- 0
for(i in 1:10000){if(test[i] == 0){a = a + 1}
}
a / 10000

[1] 0.4037

2.

用逆变换法方法生成如下概率函数的随机变量，要求写出R程序：
X 1 2 3 4 5
p 0.1 0.2 0.4 0.2 0.1

f <- function(n){{x <- rep(0,n)for (i in 1:n){u <- runif(1)if(u < 0.1) x[i] <- 1else if(u < 0.1 + 0.2) x[i] <- 2else if(u < 0.1 + 0.2 + 0.4) x[i] <- 3else if(u < 0.1 + 0.2 + 0.4 + 0.2) x[i] <- 4else x[i] <- 5}x}
}

可以测试一下：

test <- f(10000)
table(test)

test
1 2 3 4 5
1044 1980 3955 2015 1006

3.

设线性同余发生器LCG中,取M= 16,a = 5,c= 3,x0= 7.试写出一个周期{xi}的值;并问x500= ?

LCG 算法数学上基于公式：X(n + 1) = (a * X(n) + c) % m
其中, 各系数为:
模m, m > 0
系数a, 0 < a < m
增量c, 0 <= c < m
原始值(种子) 0 <= X(0) < m
其中参数c, m, a比较敏感, 或者说直接影响了伪随机数产生的质量.

据此，非常容易写出以下代码：

m <- 16
a <- 5
c <- 3
f <- function(n){{x <- rep(0,n)x[1] <- 7for (i in 1:n){x[i+1] <- (a * x[i] + c) %% m}x}
}
test <- f(500);test

测试得：

[1] 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15
[26] 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6
[51] 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9
[76] 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8
[101] 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3
[126] 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10
[151] 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13
[176] 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12
[201] 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7
[226] 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14
[251] 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1
[276] 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0
[301] 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11
[326] 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2
[351] 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5
[376] 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4
[401] 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15
[426] 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6
[451] 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9
[476] 0 3 2 13 4 7 6 1 8 11 10 5 12 15 14 9 0 3 2 13 4 7 6 1 8
[501] 11