library(dplyr)   # いろいろ
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lattice) # 3次元プロット
library(coda)    # Rhatの計算

等高線図(contour plot)を描いてみる

まずは等高線図のデータを格納するdata.frameを作成し、contour.dfと名付ける。このdata.frameにはポテンシャルエネルギーの母数(\(\theta\))thetaと運動エネルギーの母数(\(p\))pが入る。

contour.x <- seq(-10, 10, by = 0.2) 
contour.y <- seq(-10, 10, by = 0.2) 
contour.df <- expand.grid(contour.x, contour.y) # グリッド生成
colnames(contour.df) <- c("theta", "p")         # 列名の変更

グリッドが作成されたら、格納されているthetapを用いてハミルトニアンを計算する。ハミルトニアンは

\[H(\theta, p) = U(\theta) + K(p)\]

である。

ポテンシャルエネルギー

ポテンシャルエネルギ\(U(\theta)\)\(U(\theta) = mgh(\theta)\)であるが、質量(\(m\))と重力加速度(\(g\))は定数1にすると無視できるので、\(h(\theta)\)のみを用いる。この\(h(\theta)\)は、事後分布を対数変換したうえで-1を乗じたものである。ここでは\(x = 2, \sigma = 1\)の正規分布を事後分布とし、母数\(\theta\)(普通なら\(\mu\)だが…)を推定したいとする。

\[p(\theta | 2) = \frac{1}{\sqrt{2 \pi}\sigma}\exp(-\frac{(2 - \theta)^2}{2 \sigma^2})\]

\(\sigma = 1\)だし、定数も全部取りぬくとカーネルは

\[\exp(-\frac{(2 - \theta)^2}{2})\]

ここに対数と-1をかけると

\[log(\exp(-\frac{(2 - \theta)^2}{2})) \times -1 = \frac{(2 - \theta)^2}{2}\]

になる。

運動エネルギー

運動エネルギー\(K(p)\)\(K(p) = \frac{1}{2m}p^2\)であるが、ここでも質量もまた定数1とし、\(0.5 p^2\)のみを使う。

上のポテンシャルエネルギーと運動エネルギーの和がハミルトニアンである。このハミルトニアンは等高線図において「高さ」を意味することになる。

式も分かったし、早速、ハミルトニアンを計算するcontour.h関数を作成する。

# 高さを計算するcontour.h関数を作成
contour.h <- function(theta, p){
    PE <- 0.5 * (2 - theta)^2  # Potential Energy
    KE <- 0.5 * p^2            # Kinetic Energy
    return(PE + KE)            # Hamiltonian: P.E. * K.E.
}

先ほど作成したcontour.dfに新しい列hを追加し、contour.h関数を使ってハミルトニアンを計算し、格納する。

# contour.dfにh列を生成し、contour.h関数を用いて
# 等高線図における高さ(=ハミルトニアン)を格納する
contour.df <- contour.df %>% mutate(h = contour.h(theta, p))

等高線図と3次元プロットを確認する。

# 等高線図の表示
contour(x = seq(-10, 10, 0.2), y = seq(-10, 10, 0.2), 
        z = matrix(contour.df$h, nrow = 101, byrow = TRUE),
        xlab = "p", ylab = "theta")

# 3Dプロット。Latticeの使い方が分からないので
# ネットで適当にコードを持ってきました。
wireframe(h ~ theta * p, data = contour.df,
  xlab = "theta", ylab = "p", zlab = "Hamiltonian",
  main = "Surface Plot",
  drape = TRUE,
  colorkey = TRUE,
  screen = list(z = -60, x = -60)
)