#### ========== ========== ========== 分布函数与密度函数 ========== ========== ==========
rm(list = ls(all = TRUE))
graphics.off()

## ---------- 加载必要的包 ----------
library(mclust)
library(ggplot2)
library(plotly)
library(MASS)
library(fBasics)
library(LaplacesDemon)

## ---------- 钞票数据的相关矩阵 ----------
data(banknote)
round(cor(banknote[, 2:7]), digits = 2)

## ---------- X4与X6的散点图 ----------
fig = plot_ly(banknote, x = ~Bottom, y = ~Diagonal, type = "scatter", mode = "markers") %>%
  layout(plot_bgcolor = '#e5ecf6', 
         xaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, gridcolor = 'ffff'),
         yaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, gridcolor = 'ffff'))
fig

## --------- X4与X6的核密度图 ----------
fig1 = plot_ly()
dX4 = density(banknote$Bottom, kernel = "gaussian", na.rm = TRUE)
fig1 = add_lines(fig1, x = dX4$x, y = dX4$y, name = "gaussian") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, gridcolor = 'ffff'),
         yaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, gridcolor = 'ffff'))
fig1

fig2 = plot_ly()
dX6 = density(banknote$Diagonal, kernel = "gaussian")
fig2 = add_lines(fig2, x = dX6$x, y = dX6$y, name = "gaussian") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, gridcolor = 'ffff'),
         yaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, gridcolor = 'ffff'))
fig2

subplot(fig1, fig2, shareY = TRUE)

## ---------- X4与X6的密度函数的乘积的曲面图 ----------
z = outer(dX4$y, dX6$y, function(a, b) a * b)
fig3 = plot_ly(x = ~dX4$x, y = ~dX6$x, z = ~z, type = "surface")
fig3

## ---------- X4与X6的联合概率密度函数的核密度图 ----------
x = sort(banknote$Bottom)
y = sort(banknote$Diagonal)
KD_X4_6 = kde2d(x, y, n = 25)
fig4 = plot_ly(x = ~x, y = ~y, z = ~KD_X4_6$z, type = "surface")
fig4

## ---------- 多元正态分布 -----------
mu = c(3, 2)
Sigma = matrix(c(1, -1.5, -1.5, 4), nrow = 2, byrow = TRUE)
X = mvrnorm(n = 100, mu, Sigma)
colnames(X) = c("X1", "X2")
X = as.data.frame(X)

fig5 = plot_ly(X, x = ~X1, y = ~X2) %>%
  add_markers(alpha = 0.9) %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, gridcolor = 'ffff'),
         yaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, gridcolor = 'ffff'))
fig5

cor(X) # 相关矩阵

## ---------- 联合密度函数曲面 ----------
n = 80
mu1 = 3
mu2 = 2
s1 = 1
s2 = 2
rho = -1.5/2
x = seq(0, 6, length = n) * s1
y = seq(-4 , 8, length = n) * s2
f = function(x, y){
  (2 * pi * s1 * s2 * sqrt(1-rho^2))^-1 * exp(-0.5 * (1 - rho^2)^-1 * 
                   ((x-mu1)^2/s1^2 - 2 * rho * (x - mu1) * ( y - mu2) / (s1 * s2) + 
                      (y - mu2)^2/s2^2))
}
z = outer(x , y, f)

fig6 = plot_ly() %>%
  add_surface(x = ~x, y = ~y, z = ~z)
fig6

## ---------- 联合密度函数曲面的等高线 ----------
fig7 = plot_ly(x = ~x, y = ~y, z = ~z) %>%
  add_contour() 
fig7

#### ========== ========== 抽样分布与中心极限定理 ========== ==========
## ---------- Bernoulli总体容量为5、35的中心极限定理 ----------
x = array(0, dim = 1000)
for (i in 1:1000){
  y = rbinom(5, 1, prob = 1/2)
  x[i] = sqrt(5) * (mean(y) - 1/2)
}
fig8 = plot_ly()
d = density(x, kernel = 'gaussian')
fig8 = add_lines(fig8, x = d$x, y = d$y, name = "n = 5")
fig8
fig8 = add_lines(fig8, x = d$x, y = dnorm(d$x, 0, 0.5), name = "N(0, 0.25)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, gridcolor = 'ffff'),
         yaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, gridcolor = 'ffff'))
fig8

x = array(0, dim = 1000)
for (i in 1:1000){
  y = rbinom(35, 1, prob = 1/2)
  x[i] = sqrt(35) * (mean(y) - 1/2)
}
fig9 = plot_ly()
d = density(x, kernel = 'gaussian')
fig9 = add_lines(fig9, x = d$x, y = d$y, name = "n = 35")
fig9

fig9 = add_lines(fig9, x = d$x, y = dnorm(d$x, 0, 0.5), name = "N(0, 0.25)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, gridcolor = 'ffff'),
         yaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, gridcolor = 'ffff'))
fig9

subplot(fig8, fig9)

## ---------- 二维正态分布密度曲面 ----------
n = 50
mu1 = 0
mu2 = 0
s1 = 1/2
s2 = 1/2
rho = 0
x = seq(-4, 4, length = n) * s1
y = seq(-4, 4, length = n) * s2
f = function(x, y){
  (2 * pi * s1 * s2 * sqrt(1-rho^2))^-1 * exp(-0.5 * (1 - rho^2)^-1 * 
      ((x-mu1)^2/s1^2 - 2 * rho * (x - mu1) * ( y - mu2) / (s1 * s2) + (y - mu2)^2/s2^2))
}
z = outer(x , y, f)
fig10 = plot_ly(x = ~x, y = ~y, z = ~z) %>%
  add_surface() %>%
  layout(
    scene = list(
      xaxis = list(title = ""),
      yaxis = list(title = ""),
      zaxis = list(title = "")
    )
  )
fig10

## ---------- 二维核密度估计 (n=5) ----------
rm(list = ls(all = TRUE))
library(MASS)
x = array(0, dim = 1000)
y = array(0, dim = 1000)
for (i in 1:1000){
  a = rbinom(5, 1, prob = 1/2)
  x[i] = sqrt(5) * (mean(a) - 1/2)
  b = rbinom(5, 1, prob = 1/2)
  y[i] = sqrt(5) * (mean(b) - 1/2)
}
kd = kde2d(x, y, n=50, lims = c(-2, 2, -2, 2))
fig11 = plot_ly(x = ~kd$x, y = ~kd$y, z = ~kd$z) %>%
  add_surface() %>%
  layout(
    scene = list(
      xaxis = list(title = ""),
      yaxis = list(title = ""),
      zaxis = list(title = "")
    )
  )
fig11

## ---------- 二维核密度估计 (n = 85) ----------
x = array(0, dim = 1000)
y = array(0, dim = 1000)
for (i in 1:1000){
  a = rbinom(85, 1, prob = 1/2)
  x[i] = sqrt(85) * (mean(a) - 1/2)
  b = rbinom(85, 1, prob = 1/2)
  y[i] = sqrt(85) * (mean(b) - 1/2)
}
kd2 =kde2d(x, y, n=50, lims = c(-2, 2, -2, 2))
fig12 = plot_ly(x = ~kd2$x, y = ~kd2$y, z = ~kd2$z) %>%
  add_surface() %>%
  layout(
    scene = list(
      xaxis = list(title = ""),
      yaxis = list(title = ""),
      zaxis = list(title = "")
    )
  )
fig12

#### ========== ========== 厚尾分布 ========== ==========
## ---------- Cauchy 分布 ----------
x = seq(-8, 8, length = 200)
y_1 = dcauchy(x, location = 0, scale = 1)
y_2 = dcauchy(x, location = 0, scale = 1/2)
y_3 = dcauchy(x, location = 0, scale = 3)
y_4 = dcauchy(x, location = -3, scale = 2/3)
y_5 = dcauchy(x, location = 2.5, scale = 2)
fig_13 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "l = 0, s = 1") %>%
  add_lines(x = ~x, y = ~y_2, name = "l = 0, s = 1/2") %>%
  add_lines(x = ~x, y = ~y_3, name = "l = 0, s = 3") %>%
  add_lines(x = ~x, y = ~y_4, name = "l = -1, s = 2/3") %>%
  add_lines(x = ~x, y = ~y_5, name = "l = 2.5, s = 2")
fig_13 = fig_13 %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = '#ffff', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""))
fig_13

## ---------- Cauchy 分布是厚尾分布 ----------
x = seq(-6, 6, length = 200)
y_1 = dnorm(x, mean = 0, sd = 1)
y_2 = dcauchy(x, location = 0, scale = 1)
fig_14 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "Normal") %>%
  add_lines(x = ~x, y = ~y_2, name = "Cauchy")
fig_14 = fig_14 %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""))
fig_14

## ---------- 双曲线 (HYP) 分布 ----------
# ----- 双曲线分布的概率密度曲线 -----
x = seq(-6, 6, length = 200)
y_1 = dhyp(x, alpha = 1, beta = 0, delta = 1, mu = 0)
y_2 = dhyp(x, alpha = 1, beta = 0, delta = 1/2, mu = 0)
y_3 = dhyp(x, alpha = 1, beta = 0, delta = 1/100, mu = 0)
y_4 = dhyp(x, alpha = 1, beta = 0, delta = 4, mu = 0)
fig_15 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "HYP (delta = 1)") %>%
  add_lines(x = ~x, y = ~y_2, name = "HYP (delta = 1/2)") %>%
  add_lines(x = ~x, y = ~y_3, name = "HYP (delta = 1/100)") %>%
  add_lines(x = ~x, y = ~y_4, name = "HYP (delta = 4)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "双曲线分布的概率密度函数"))
fig_15

# ----- 双曲线分布的分布函数曲线 -----
x = seq(-6, 6, length = 200)
y_1 = phyp(x, alpha = 1, beta = 0, delta = 1, mu = 0)
y_2 = phyp(x, alpha = 1, beta = 0, delta = 1/2, mu = 0)
y_3 = phyp(x, alpha = 1, beta = 0, delta = 1/100, mu = 0)
y_4 = phyp(x, alpha = 1, beta = 0, delta = 4, mu = 0)
fig_16 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "HYP (delta = 1)") %>%
  add_lines(x = ~x, y = ~y_2, name = "HYP (delta = 1/2)") %>%
  add_lines(x = ~x, y = ~y_3, name = "HYP (delta = 1/100)") %>%
  add_lines(x = ~x, y = ~y_4, name = "HYP (delta = 4)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "双曲线分布的概率分布函数"))
fig_16

## ---------- 正态逆高斯 (NIG) 分布 ----------
# ----- 正态逆高斯分布的概率密度曲线 -----
x = seq(-4, 4, length = 200)
y_1 = dnig(x, alpha = 1, beta = 0, delta = 1, mu = 0)
y_2 = dnig(x, alpha = 1, beta = 0, delta = 1/2, mu = 0)
y_3 = dnig(x, alpha = 1, beta = 0, delta = 1/5, mu = 0)
y_4 = dnig(x, alpha = 1, beta = 0, delta = 4, mu = 0)
fig_17 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "NIG (delta = 1)") %>%
  add_lines(x = ~x, y = ~y_2, name = "NIG (delta = 1/2)") %>%
  add_lines(x = ~x, y = ~y_3, name = "NIG (delta = 1/5)") %>%
  add_lines(x = ~x, y = ~y_4, name = "NIG (delta = 4)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "正态逆高斯分布的概率密度函数"))
fig_17

# ----- 正态逆高斯分布的分布函数曲线 -----
x = seq(-4, 4, length = 200)
y_1 = pnig(x, alpha = 1, beta = 0, delta = 1, mu = 0)
y_2 = pnig(x, alpha = 1, beta = 0, delta = 1/2, mu = 0)
y_3 = pnig(x, alpha = 1, beta = 0, delta = 1/5, mu = 0)
y_4 = pnig(x, alpha = 1, beta = 0, delta = 4, mu = 0)
fig_18 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "NIG (delta = 1)") %>%
  add_lines(x = ~x, y = ~y_2, name = "NIG (delta = 1/2)") %>%
  add_lines(x = ~x, y = ~y_3, name = "NIG (delta = 1/5)") %>%
  add_lines(x = ~x, y = ~y_4, name = "NIG (delta = 4)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "正态逆高斯分布的概率分布函数"))
fig_18

## ---------- 广义双曲线 (GH) 分布 ----------
# ----- 广义双曲线分布的概率密度函数 -----
x = seq(-4, 4, length = 200)
y_1 = dgh(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = -1/2)
y_2 = dgh(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1/2)
y_3 = dgh(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 2)
fig_19 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "GH (lambda = -1/2)") %>%
  add_lines(x = ~x, y = ~y_2, name = "GH (lambda = 1/2)") %>%
  add_lines(x = ~x, y = ~y_3, name = "GH (lambda = 2)") %>%
    layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "广义双曲线分布的概率密度函数"))
fig_19

# ----- 广义双曲线分布的概率分布函数 -----
x = seq(-4, 4, length = 200)
y_1 = pgh(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = -1/2)
y_2 = pgh(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1/2)
y_3 = pgh(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 2)
fig_20 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "GH (lambda = -1/2)") %>%
  add_lines(x = ~x, y = ~y_2, name = "GH (lambda = 1/2)") %>%
  add_lines(x = ~x, y = ~y_3, name = "GH (lambda = 2)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "广义双曲线分布的概率分布函数"))
fig_20

# ----- 不同概率密度函数的对比 -----
x = seq(-6, 6, length = 200)
y_1 = dgh(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 0.5)
y_2 = dhyp(x, alpha = 1, beta = 0, delta = 1, mu = 0)
y_3 = dnig(x, alpha = 1, beta = 0, delta = 1, mu = 0)
fig_21 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "GH (lambda = 1/2)") %>%
  add_lines(x = ~x, y = ~y_2, name = "HYP") %>%
  add_lines(x = ~x, y = ~y_3, name = "NIG") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "概率密度函数"))
fig_21

# ----- 不同概率分布函数的对比 -----
x = seq(-6, 6, length = 200)
y_1 = pgh(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 0.5)
y_2 = phyp(x, alpha = 1, beta = 0, delta = 1, mu = 0)
y_3 = pnig(x, alpha = 1, beta = 0, delta = 1, mu = 0)
fig_22 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "GH (lambda = 1/2)") %>%
  add_lines(x = ~x, y = ~y_2, name = "HYP") %>%
  add_lines(x = ~x, y = ~y_3, name = "NIG") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "概率分布函数"))
fig_22

## ---------- 学生 t 分布 ----------
# ----- t 分布的密度函数曲线 -----
x = seq(-6, 6, length = 200)
y_1 = dt(x, 15)
y_2 = dt(x, 5)
y_3 = dt(x, 1)
y_4 = dnorm(x, 0, 1)
fig_23 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "t (n = 15)") %>%
  add_lines(x = ~x, y = ~y_2, name = "t (n = 5)") %>%
  add_lines(x = ~x, y = ~y_3, name = "t (n = 1)") %>%
  add_lines(x = ~x, y = ~y_4, name = "N(0, 1)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "概率密度曲线"))
fig_23

# ----- t 分布的分布函数曲线 -----
x = seq(-6, 6, length = 200)
y_1 = pt(x, 15)
y_2 = pt(x, 5)
y_3 = pt(x, 1)
y_4 = pnorm(x, 0, 1)
fig_24 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "t (n = 15)") %>%
  add_lines(x = ~x, y = ~y_2, name = "t (n = 5)") %>%
  add_lines(x = ~x, y = ~y_3, name = "t (n = 1)") %>%
  add_lines(x = ~x, y = ~y_4, name = "N(0, 1)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "分布函数曲线"))
fig_24

# ----- t 分布尾部的比较 -----
x = seq(2.5, 5, length = 100)
y_1 = dt(x, 1)
y_2 = dt(x, 3)
y_3 = dt(x, 9)
y_4 = dt(x, 45)
y_5 = dnorm(x, 0, 1)
fig_25 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "t (n = 1)") %>%
  add_lines(x = ~x, y = ~y_2, name = "t (n = 3)") %>%
  add_lines(x = ~x, y = ~y_3, name = "t (n = 9)") %>%
  add_lines(x = ~x, y = ~y_4, name = "t (n = 45)") %>%
  add_lines(x = ~x, y = ~y_5, name = "N(0, 1)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "t 分布密度曲线尾部的比较"))
fig_25

## ---------- Laplace 分布 ----------
# ----- Laplace 分布的密度函数曲线 -----
x = seq(-6, 6, length = 201)
y_1 = dlaplace(x, location = 0, scale = 1)
y_2 = dlaplace(x, location = 0, scale = 1.5)
y_3 = dlaplace(x, location = 0, scale = 2)
fig_26 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "L (theta = 1)") %>%
  add_lines(x = ~x, y = ~y_2, name = "L (theta = 1.5)") %>%
  add_lines(x = ~x, y = ~y_3, name = "L (theta = 2)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "Laplace 分布密度曲线"))
fig_26

# ----- Laplace 分布的分布函数曲线 -----
x = seq(-6, 6, length = 201)
y_1 = plaplace(x, location = 0, scale = 1)
y_2 = plaplace(x, location = 0, scale = 1.5)
y_3 = plaplace(x, location = 0, scale = 2)
fig_27 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "L (theta = 1)") %>%
  add_lines(x = ~x, y = ~y_2, name = "L (theta = 1.5)") %>%
  add_lines(x = ~x, y = ~y_3, name = "L (theta = 2)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "Laplace 分布的分布函数曲线"))
fig_27

## ---------- Cauchy 分布 ----------
# ----- Cauchy 分布的密度函数曲线 -----
x = seq(-6, 6, length = 201)
y_1 = dcauchy(x, location = 0, scale = 1)
y_2 = dcauchy(x, location = 0, scale = 1.5)
y_3 = dcauchy(x, location = 0, scale = 2)
fig_28 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "Cauchy (s = 1)") %>%
  add_lines(x = ~x, y = ~y_2, name = "Cauchy (s = 1.5)") %>%
  add_lines(x = ~x, y = ~y_3, name = "Cauchy (s = 2)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "Cauchy 分布的密度函数曲线"))
fig_28

# ----- Cauchy 分布的分布函数曲线 -----
x = seq(-6, 6, length = 201)
y_1 = pcauchy(x, location = 0, scale = 1)
y_2 = pcauchy(x, location = 0, scale = 1.5)
y_3 = pcauchy(x, location = 0, scale = 2)
fig_29 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "Cauchy (s = 1)") %>%
  add_lines(x = ~x, y = ~y_2, name = "Cauchy (s = 1.5)") %>%
  add_lines(x = ~x, y = ~y_3, name = "Cauchy (s = 2)") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "Cauchy 分布的分布函数曲线"))
fig_29


#### ---------- ---------- 混合模型 ---------- ----------
## ----- 80% 的 N(0, 1) 与 20% 的 N(0, 9) 构成的正态混合分布的密度函数曲线 -----
x = seq(-6, 6, length = 201)
y_1 = dnorm(x, mean = 0, sd = 1)
y_2 = dnorm(x, mean = 0, sd = 3)
y_3 = 0.8 * y_1 + 0.2 * y_2
fig_30 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "N(0,1)") %>%
  add_lines(x = ~x, y = ~y_2, name = "N(0,9)") %>%
  add_lines(x = ~x, y = ~y_3, name = "GM") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "正态混合模型的密度函数曲线"))
fig_30

# ----- 80% 的 N(0, 1) 与 20% 的 N(0, 9) 构成的正态混合分布的分布函数曲线 -----
x = seq(-6, 6, length = 201)
y_1 = pnorm(x, mean = 0, sd = 1)
y_2 = pnorm(x, mean = 0, sd = 3)
y_3 = 0.8 * y_1 + 0.2 * y_2
fig_31 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "N(0,1)") %>%
  add_lines(x = ~x, y = ~y_2, name = "N(0,9)") %>%
  add_lines(x = ~x, y = ~y_3, name = "GM") %>%
  layout(plot_bgcolor='#e5ecf6',
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "正态混合模型的分布函数曲线"))
fig_31

## ---------- 广义双曲分布 ----------
# ----- 广义双曲分布以指数速度衰减 -----
x = seq(4, 6, length = 101)
y_1 = dgh(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 0.5)
y_2 = dgh(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1.5)
y_3 = dnig(x, alpha = 1, beta = 0, delta = 1, mu = 0)
y_4 = dhyp(x, alpha = 1, beta = 0, delta = 1, mu = 0)
fig_32 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = TeX("\\text{GH} (\\lambda = 0.5)")) %>%
  add_lines(x = ~x, y = ~y_2, name = TeX("\\text{GH} (\\lambda = 1.5)")) %>%
  add_lines(x = ~x, y = ~y_3, name = "NIG") %>%
  add_lines(x = ~x, y = ~y_4, name = "HYP") %>%
  layout(plot_bgcolor='#e5ecf6', title = "Tail comparison of GH distributionns (pdf)",
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = "GH 的密度曲线")) %>%
  config(mathjax = "cdn")
fig_32

# ----- 正态逆高斯(NIG)、Laplace、Cauchy、标准正态(Gaussian) 分布的密度曲线的比较 -----
x = seq(-6, 6, length = 201)
y_1 = dnig(x, alpha = 1, beta = 0, delta = 1, mu = 0)
y_2 = dlaplace(x, location = 0, scale = 1)
y_3 = dnorm(x, 0, 1)
y_4 = dcauchy(x, location = 0, scale = 1)
fig_33 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "NIG") %>%
  add_lines(x = ~x, y = ~y_2, name = "Laplace") %>%
  add_lines(x = ~x, y = ~y_3, name = "Gaussian") %>%
  add_lines(x = ~x, y = ~y_4, name = "Cauchy") %>%
  layout(plot_bgcolor='#e5ecf6', title = "Distributions Comparison",
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""))
fig_33

# ----- 正态逆高斯(NIG)、Laplace、Cauchy、标准正态(Gaussian) 分布尾部的比较 -----
x = seq(-6, -4, length = 101)
y_1 = dnig(x, alpha = 1, beta = 0, delta = 1, mu = 0)
y_2 = dlaplace(x, location = 0, scale = 1)
y_3 = dnorm(x, 0, 1)
y_4 = dcauchy(x, location = 0, scale = 1)
fig_34 = plot_ly() %>%
  add_lines(x = ~x, y = ~y_1, name = "NIG") %>%
  add_lines(x = ~x, y = ~y_2, name = "Laplace") %>%
  add_lines(x = ~x, y = ~y_3, name = "Gaussian") %>%
  add_lines(x = ~x, y = ~y_4, name = "Cauchy") %>%
  layout(plot_bgcolor='#e5ecf6', title = "Tail Comparison",
         xaxis = list(zerolinecolor = 'cyan', zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""),
         yaxis = list(zerolinecolor = "cyan", zerolinewidth = 2, 
                      gridcolor = 'ffff', title = ""))
fig_34

#### ========== ========== 自助法 Bootstrap ========== ==========

## ----- 经验分布函数(n=10) -----
n = 10
x = rnorm(n, 0, 1)
plot(ecdf(x), xlab = 'X', ylab ='', xlim = c(-3, 3), main = 'ECDF and CDF', 
     verticals = TRUE, do.points = FALSE)
curve(pnorm(x, 0, 1), -3, 3, add = TRUE, col = 'red', lwd = 2)

# 经验分布函数(n=100)
n = 100
x = rnorm(n, 0, 1)
plot(ecdf(x), xlab = 'X', ylab ='', xlim = c(-3, 3), main = 'ECDF and CDF', 
     verticals = TRUE, do.points = FALSE)
curve(pnorm(x, 0, 1), -3, 3, add = TRUE, col = 'red', lwd = 2)

# 经验分布函数(n=1000)
n = 1000
x <- rnorm(n, 0, 1)
plot(ecdf(x), xlab = 'X', ylab ='', xlim = c(-3, 3), main = 'ECDF and CDF', 
     verticals = TRUE, do.points = FALSE)
curve(pnorm(x, 0, 1), -3, 3, add = TRUE, col = 'red', lwd = 2)

# 自助样本(n=100)
n = 100
x = rnorm(n, 0, 1)
y1 = sample(x, n, replace = TRUE)
y2 = sample(x, n, replace = TRUE)
plot(ecdf(x), xlab = 'X', ylab ='', xlim = c(-3, 3), 
     main = 'ECDF and 2 bootstrap ECDFs, n = 100',
     verticals = TRUE, do.points = FALSE, lwd = 2, col = 'black')
plot(ecdf(y1), verticals = TRUE, do.points = FALSE, col = 'red', add = TRUE)
plot(ecdf(y2), verticals = TRUE, do.points = FALSE, col = 'blue', add = TRUE)
curve(pnorm(x, 0, 1), -3, 3, add = TRUE, col = 'cyan3', lwd = 2)