6-9，最好能用R做

Question

是大逆不道不好看们的奶奶手机奥斯卡上看完子子孙孙本山那是宝宝睡吧

Answer 1

引自强大而免费的微信小程序：皆我百晓生。

6.使用R语言计算P{X>30}，其中X~B(40,0.6)：

# 使用二项分布函数计算
p <- 1 - pbinom(30, size=40, prob=0.6)
p

通过二项分布函数计算得到的概率精度较高。

使用正态近似方法计算P{X>30}：

mu <- 40 * 0.6
sigma <- sqrt(40 * 0.6 * 0.4)
p <- 1 - pnorm(30, mean=mu, sd=sigma)
p

通过正态近似方法计算得到的概率会有一定的误差。

7.生成两个随机数发生器：

# 混合同余法1
hx <- function() {
  m <- 2^31 - 1
  a <- 16807
  b <- 0
  hx <- function(x) {
    x <- (a * x + b) %% m
    return(x)
  }
  return(hx)
}

# 混合同余法2
hy <- function() {
  m <- 2^31 - 1
  a <- 40014
  b <- 0
  hy <- function(y) {
    y <- (a * y + b) %% m
    return(y)
  }
  return(hy)
}

# 组合同余发生器
hz <- function(x, y) {
  m <- 2^31 - 1
  hz <- function() {
    x <- hx(x)
    y <- hy(y)
    z <- (x - y) %% m
    return(z/m)
  }
  return(hz)
}

# 测试生成随机数
seed_x <- 1234
seed_y <- 5678

generator <- hz(seed_x, seed_y)

# 生成10个随机数
random_numbers <- replicate(10, generator(), simplify = TRUE)
random_numbers

8.画出密度曲线和分布曲线，并计算分位数和生成随机数：

# 密度曲线和分布曲线
x <- seq(-15, 21, length.out=100)
density <- dnorm(x, mean=3, sd=4)
cumulative <- pnorm(x, mean=3, sd=4)

plot(x, density, type="l", main="Normal Distribution", xlab="x", ylab="Density")
lines(x, cumulative, col="red")
legend("topleft", c("Density", "Cumulative"), lty=1, col=c("black", "red"))

# 分位数
p25 <- qnorm(0.25, mean=3, sd=4)
p50 <- qnorm(0.5, mean=3, sd=4)
p75 <- qnorm(0.75, mean=3, sd=4)
p25
p50
p75

# 生成40个随机数
random_numbers <- rnorm(40, mean=3, sd=4)
random_numbers

9.使用逆变换法、舍选法和复合法生成随机数并比较效率：

# 逆变换法生成随机数
inverse_transform <- function(n) {
  u <- runif(n)
  x <- sqrt(2 * u)
  return(x)
}

# 舍选法生成随机数
rejection_sampling <- function(n) {
  x <- NULL
  while(length(x) < n) {
    u <- runif(1)
    y <- runif(1)
    if (y <= 0.5 + x) {
      x <- c(x, u)
    }
  }
  return(x)
}

# 复合法生成随机数
composite_method <- function(n) {
  x <- NULL
  while(length(x) < n) {
    u <- runif(1)
    y <- runif(1)
    if (y <= 1/2 + u^2/2) {
      x <- c(x, u)
    }
  }
  return(x)
}

# 生成1000个随机数并比较效率
library(microbenchmark)

samples_inv <- microbenchmark(inverse_transform(1000))
samples_rej <- microbenchmark(rejection_sampling(1000))
samples_com <- microbenchmark(composite_method(1000))

names(samples_inv) <- "Inverse Transform"
names(samples_rej) <- "Rejection Sampling"
names(samples_com) <- "Composite Method"

summary(samples_inv)
summary(samples_rej)
summary(samples_com)