“应用统计田恬 统计计算作业一”

第 1 题

（1）

令 \(u=F(x)\)，解\(x\)可得

\[ F^{-1}(u)=\frac{b}{(1-u)^{1/a}} =b(1-u)^{-1/a}, \qquad 0≤u<1. \]

（2）

当 \(a=2\)、\(b=2\) 时，逆变换函数为

\[ F^{-1}(u)=\frac{2}{\sqrt{1-u}}. \]

set.seed(111)

n <- 1000
a <- 2
b <- 2

u <- runif(n)
x <- b / (1 - u)^(1 / a)

（3）

Pareto\((a,b)\) 的密度函数为

\[ f(x)=\frac{ab^a}{x^{a+1}}, \qquad x\ge b. \]

所以 Pareto\((2,2)\) 的密度为 \(f(x)=8/x^3\)（\(x\ge 2\)）。

h <- hist(
  x,
  prob = TRUE,
  breaks = "FD",
  xlim = c(2, 30),
  ylim = c(0, 1.05),
  main = bquote(f(x) == a*b^a/x^{a+1})
)
y <- seq(2, 30, 0.01)
lines(y, 8/y^3)

第 2 题

（1）

set.seed(222)

n <- 1000
u <- runif(n)

x_inverse <- integer(n)
x_inverse[0.0 <= u & u < 0.1] <- 0
x_inverse[0.1 <= u & u < 0.3] <- 1
x_inverse[0.3 <= u & u < 0.5] <- 2
x_inverse[0.5 <= u & u < 0.7] <- 3
x_inverse[0.7 <= u & u <= 1.0] <- 4

（2）

x_value <- 0:4
prob <- c(0.1, 0.2, 0.2, 0.2, 0.3)

x <- sample(
  x_value,
  size = n,
  replace = TRUE,
  prob = prob
)

（3）

round(rbind(table(x)/n,prob),3)

##          0     1     2     3     4
##      0.099 0.212 0.189 0.184 0.316
## prob 0.100 0.200 0.200 0.200 0.300

第 3 题

（1）

Beta(3,2) 的密度函数为 \(12x^2(1-x)\)，令 \(g(x)\) 为 U(0,1) 分布的密度函数。这样对于 \(0<x<1\) 都有 \(f(x)/g(x)\le 12\)， 因此可以令 \(c=12\)。

n=1000
k=0
j=0
y=numeric(n)
while(k<n){
  u=runif(1)
  j=j+1
  x=runif(1)
  if (x^2*(1-x)>u){
    k=k+1
    y[k]=x
  }
}
j

## [1] 11859

（2）

hist(
  y,
  prob = TRUE,
  breaks = "FD",
  xlim = c(0, 1),
  ylim = c(0, 2),
  main = bquote(f(x) == 12*x^2*(1-x))
)
x <- seq(0, 1, by = 0.01)
y <- 12*x^2*(1-x)
lines(x, y, lwd = 2)

第 4 题

（1）

lognorm<-function(n,mu,sigma2){
  u<-runif(n)
  z<-qnorm(u)
  x<-exp(mu+sqrt(sigma2)*z)
  return(x)
}

（2）

set.seed(111)
y=lognorm(1000,0,1)

（3）

x_right=unname(quantile(y,0.99))
hist(
  y,
  prob = TRUE,
  breaks = "FD",
  xlim = c(0, x_right),
  ylim = c(0, 0.8)
)
z1=seq(0,x_right,length.out=1000)
z2=dlnorm(z1,0,1)
lines(z1,z2,lwd=2)

第 5 题

loc.mix<-function(n,p,mu1,mu2,sigma2){
  n1<-rbinom(1,size=n,prob=p)
  n2<-n-n1
  x1<-rnorm(n1,mu1,sqrt(sigma2))
  x2<-rnorm(n2,mu2,sqrt(sigma2))
  X<-c(x1,x2)
  return(sample(X))
}
set.seed(111)
x1=loc.mix(1000,.1,0,3,1)
x2=loc.mix(1000,.3,0,3,1)
x3=loc.mix(1000,.5,0,3,1)
x4=loc.mix(1000,.7,0,3,1)
x5=loc.mix(1000,.9,0,3,1)

第 6 题

rmvn.eigen<-function(n,mu,Sigma){
  d<-length(mu)
  ev<-eigen(Sigma,symmetric=TRUE)
  lambda<-ev$values
  V<-ev$vectors
  R<-V%*%diag(sqrt(lambda)) %*% t(V)
  Z<-matrix(rnorm(n*d),nrow=n,ncol=d)
  X<-Z%*%R+matrix(mu,n,d,byrow=TRUE)
  X
}
rmvn.svd<-function(n,mu,Sigma){
  d<-length(mu)
  S<-svd(Sigma)
  R<-S$u%*%diag(sqrt(S$d)) %*% t(S$v)
  Z<-matrix(rnorm(n*d),nrow=n,ncol=d)
  X<-Z%*%R+matrix(mu,n,d,byrow=TRUE)
  X  
}
rmvn.Choleski<-function(n,mu,Sigma){
  d<-length(mu)
  Q<-chol(Sigma)
  Z<-matrix(rnorm(n*d),nrow=n,ncol=d)
  X<-Z%*%Q+matrix(mu,n,d,byrow=TRUE)
  X  
}
set.seed(111)

mu <- c(0, 1, 2)
Sigma <- matrix(
  c(1.0, -0.5,  0.5,
   -0.5,  1.0, -0.5,
    0.5, -0.5,  1.0),
  nrow = 3,
  ncol = 3,
  byrow = TRUE
)

X=rmvn.eigen(200,mu,Sigma)
pairs(X)

Y=rmvn.svd(200,mu,Sigma)
pairs(Y)

Z=rmvn.Choleski(200,mu,Sigma)
pairs(Z)