# https://w...content-available-to-author-only...t.cc/bbs/R_Language/M.1540859231.A.848.html

theta <- 8 # 本例以 theta = 8 為例
n.delta <- 500 # 設定方格在 x 和 y 的格數
x.sim <- seq(-theta, theta, length = n.delta)
y.sim <- seq(-14, 14, length = n.delta) # 太窄或太寬可以調整
f.xy <-
  function(x, y, th) {
    1 / (2 * th * sqrt(2 * pi)) * exp(-1 / 2 * (y - x) ^ 2)
  }
# 求 Pxy(x,y): z.sim 為 200 * 200 方格
z.sim <- outer(x.sim, y.sim, FUN = f.xy, th = theta)

# 求 Py(y)，已知 Px 是均勻分配故可簡化算法
Py <- apply(z.sim, 2, sum) / n.delta * (2 * theta)

# trapezoidal integration 檢查 Py dy 積分
sum((Py[1:(length(Py) - 1)] + Py[2:(length(Py))]) * diff(y.sim)[1] / 2)

# Pxy 視覺化
windows(5, 5)
surf3D(
  matrix(rep(x.sim, n.delta), ncol = n.delta),
  matrix(rep(y.sim, each = n.delta), ncol = n.delta),
  z.sim,
  bty = "f",
  theta = 110,
  phi = 30,
  r = 3,
  ticktype = "detailed"
)

# Py 視覺化
windows(5, 5)
plot(y.sim,
     Py,
     xlab = "y",
     ylab = "P_y(y)",
     type = "l")