# 第 11 章 三维可视化

library(barsurf)
library(plotrix)
library(scatterplot3d)
library(plot3D)
library(MBA)

## 11.5 透视图

x <- seq(-10, 10, length = 30)
y <- x
f <- function(x, y) {
r <- sqrt(x^2 + y^2)
10 * sin(r) / r
}
z <- outer(x, y, f)
z[is.na(z)] <- 1
op <- par(bg = "white")
nrz <- nrow(z)
ncz <- ncol(z)
jet.colors <- colorRampPalette(c("gray80", "gray10"))
nbcol <- 100
color <- jet.colors(nbcol)
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
facetcol <- cut(zfacet, nbcol)

persp(x, y, z,
theta = 30, phi = 30,
expand = 0.5, col = color[facetcol]
)

persp(x, y, z, xaxs = "i", expand = 0.5, phi = 20, theta = 60, col = color[facetcol])

persp(x, y, z,
theta = 45, phi = 20,
expand = 0.5, col = color[facetcol],
r = 180,
ltheta = 120,
# ticktype = "detailed", # 坐标轴上刻度数字
# box=FALSE, # 长方体框线
# nticks=6, # 刻度间隔数目
xlab = "X", ylab = "Y", zlab = "Sinc( r )"
# border=30
)

theta参数给出了主要方向，控制三维图的左右，phi给出纬度，expand 控制三维图的立体性

x <- seq(-3, 3, by = 0.25)
y <- seq(-3, 3, by = 0.25)
d <- expand.grid(x = x, y = y)
z <- c(data = NA, 1089)
b0 <- 5.628
b1 <- 0
b2 <- 0
b3 <- -.1
b4 <- .1
b5 <- -.1
k <- 1
for (i in 1:25) {
for (j in 1:25) {
z[k] <- b0 + b1 * x[i] + b2 * y[j] + b3 * x[i] * x[i] + b4 * x[i] * y[j] + b5 * y[j] * y[j]
k <- k + 1
}
}
library(rsm)
data.lm <- lm(z ~ poly(x, y, degree = 2), data = d)
persp(data.lm, x ~ y,
zlim = c(0, max(z)),
contour = list(z = "bottom", col = "colors"), theta = -55, phi = 25
)

res1 <- persp(data.lm, x ~ y,
zlim = c(0, max(z)),
contour = list(z = "bottom", col = "colors"), theta = -55, phi = 25
)
xy <- matrix(c((-3 - 8) / 5, -3, (3 - 8) / 5, 3), ncol = 2, byrow = T)
lines(trans3d(xy[, 2], xy[, 1], 0, pmat = res1$y ~ x$transf), col = 3)

set.seed(1234)
n <- 20 # 随机数的个数
x <- rexp(n, rate = 5)
m <- 40 # 网格数
mv <- seq(mean(x) - 1.5 * sd(x) / sqrt(n),
mean(x) + 1.5 * sd(x) / sqrt(n),
length.out = m
) # mu 均值范围
sv <- seq(0.8 * sd(x), 1.5 * sd(x), length.out = m) # 标准差的范围
z <- matrix(NA, m, m)
loglikelihood <- function(b) {
-sum(dnorm(x, b[1], b[2], log = TRUE))
}

for (i in 1:m) {
for (j in 1:m) {
z[i, j] <- -loglikelihood(c(mv[i], sv[j]))
}
}

nbcol <- 100
color <- hcl.colors(nbcol)
zfacet <- z[-1, -1] + z[-1, -m] + z[-m, -1] + z[-m, -m]
facetcol <- cut(zfacet, nbcol)
# "\n" adds one line before the label
persp(mv, sv, z,
xlab = "\n mu", ylab = "\n sigma", zlab = "\n log-likelihood",
phi = 35, theta = -30, col = color[facetcol]
)

library(latex2exp)
# 代码来自  http://www.ejwagenmakers.com/misc/Plotting_3d_in_R.pdf
mu1 <- 0 # setting the expected value of x1
mu2 <- 0 # setting the expected value of x2
s11 <- 10 # setting the variance of x1
s12 <- 15 # setting the covariance between x1 and x2
s22 <- 10 # setting the variance of x2
rho <- 0.5 # setting the correlation coefficient between x1 and x2
x1 <- seq(-10, 10, length = 41) # generating the vector series x1
x2 <- x1 # copying x1 to x2
# setting up the function of the multivariate normal density
f <- function(x1, x2) {
term1 <- 1 / (2 * pi * sqrt(s11 * s22 * (1 - rho^2)))
term2 <- -1 / (2 * (1 - rho^2))
term3 <- (x1 - mu1)^2 / s11
term4 <- (x2 - mu2)^2 / s22
term5 <- -2 * rho * ((x1 - mu1) * (x2 - mu2)) / (sqrt(s11) * sqrt(s22))
term1 * exp(term2 * (term3 + term4 - term5))
}
z <- outer(x1, x2, f) # calculating the density values
nrz <- nrow(z)
ncz <- ncol(z)
nbcol <- 100
color <- hcl.colors(100)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)

par(mar = c(4.1, 4.1, 4.5, 1.5), ps = 10)
persp(x1, x2, z,
xlab = "\n x1",
ylab = "\n x2",
zlab = "\n\n f(x1,x2)",
# xlab = TeX('$x_{1}$'), # latex2exp 其实是使用 LaTeX 语法将 LaTeX 公式翻译为 R 能接受的表达式形式
# ylab = TeX('$x_{2}$'),
# zlab = TeX('$f(x_{1},x_{2})$'),
main = "Two dimensional Normal Distribution",
col = color[facetcol], border = NA, theta = 30, phi = 20,
r = 50, d = 0.1, expand = 0.5, ltheta = 90, lphi = 180,
shade = 0.1, ticktype = "detailed", nticks = 5, box = TRUE
)
mtext(expression(list(
mu[1] == 0, mu[2] == 0, sigma[11] == 10,
sigma[22] == 10, sigma[12] == 15, rho == 0.5
)),
side = 3
)

mtext(expression(italic(f) ~ group("(", list(x[1], x[2]), ")") == frac(1, 2 ~ pi ~ sqrt(sigma[11] ~ sigma[22] ~ (1 - rho^2))) ~ exp ~
bgroup(
"{",
paste(
-frac(1, 2(1 - rho^2)) * phantom(0),
bgroup(
"[",
frac((x[1] ~ -~ mu[1])^2, sigma[11]) ~
-~2 ~ rho ~ frac(x[1] ~ -~ mu[1], sqrt(sigma[11])) ~ frac(x[2] ~ -~ mu[2], sqrt(sigma[22])) ~
+~ frac((x[2] ~ -~ mu[2])^2, sigma[22]),
"]"
)
),
"}"
)), side = 1, line = 3)

# 代码来自  http://www.ejwagenmakers.com/misc/Plotting_3d_in_R.pdf
mu1 <- 0 # setting the expected value of x1
mu2 <- 0 # setting the expected value of x2
s11 <- 10 # setting the variance of x1
s12 <- 15 # setting the covariance between x1 and x2
s22 <- 10 # setting the variance of x2
rho <- 0.5 # setting the correlation coefficient between x1 and x2
x1 <- seq(-10, 10, length = 41) # generating the vector series x1
x2 <- x1 # copying x1 to x2
f <- function(x1, x2) {
term1 <- 1 / (2 * pi * sqrt(s11 * s22 * (1 - rho^2)))
term2 <- -1 / (2 * (1 - rho^2))
term3 <- (x1 - mu1)^2 / s11
term4 <- (x2 - mu2)^2 / s22
term5 <- -2 * rho * ((x1 - mu1) * (x2 - mu2)) / (sqrt(s11) * sqrt(s22))
term1 * exp(term2 * (term3 + term4 - term5))
} # setting up the function of the multivariate normal density
z <- outer(x1, x2, f) # calculating the density values
nrz <- nrow(z)
ncz <- ncol(z)
nbcol <- 100
color <- hcl.colors(100)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)
par(mar = c(4.1, 4.1, 4.5, 1.5))
persp(x1, x2, z,
xlab = "$x_{1}$",
ylab = "$x_{2}$",
zlab = "$f(x_{1},x_{2})$",
main = "Two dimensional Normal Distribution",
col = color[facetcol], border = NA, theta = 30, phi = 20,
r = 50, d = 0.1, expand = 0.5, ltheta = 90, lphi = 180,
shade = 0.1, ticktype = "detailed", nticks = 5, box = TRUE
)
mtext("$\\mu_1 = 0,\\mu_2 = 0,\\sigma_{11} = 10,\\sigma_{22} = 10,\\sigma_{12} = 15, \\rho = 0.5$", side = 3)
mtext("$f(x_{1},x_{2}) = \\frac{1}{2\\pi\\sqrt{\\sigma_{11}\\sigma_{22}(1-\\rho^2)}}\\exp\\big\\{-\\frac{1}{2(1-\\rho^2)}[\\frac{(x_1 - \\mu_1)^2}{\\sigma_{11}} - 2\\rho\\frac{(x_1 - \\mu_1)(x_2 - \\mu_2)}{\\sqrt{\\sigma_{11}}\\sqrt{\\sigma_{22}}} + \\frac{(x_2 - \\mu_2)^2}{\\sigma_{22}}]\\big\\}$",
side = 1, line = 2, cex = 1.5
)
library(lattice)
wireframe(z ~ x1 + x2,
data = data.frame(x1 = x1, x2 = rep(x2, each = length(x1)), z = z),
xlab = expression(x[1]), ylab = expression(x[2]),
zlab = expression(italic(f) ~ group("(", list(x[1], x[2]), ")")),
colorkey = TRUE, drape = TRUE
)
## volcano  ## 87 x 61 matrix
wireframe(volcano,
aspect = c(61 / 87, 0.4),
light.source = c(10, 0, 10)
)

# https://stackoverflow.com/questions/41190525/adjust-margins-in-persp-persp3d-in-r
# https://stackoverflow.com/questions/37571376/how-to-customize-nticks-in-persp-r
# https://stackoverflow.com/questions/43507680/are-xlab-ylab-and-zlab-in-persp-incompatible-with-bquote
x <- seq(-10, 10, len = 30)
y <- seq(0, 5, len = 30)
f <- function(x, y) {
dnorm(2, x, y)
}
z <- outer(x, y, f)
persp(x, y, z,
theta = 30, phi = 30, expand = 0.5, col = "lightblue",
xlab = "\u03bc", ylab = "\u03c3\u207F",
zlab = paste("Likelihood ", "(\u03bc,\u03c3\u207F)", sep = "")
)

library(MBA)
data(LIDAR) # 一小部分光探测和测距雷达数据，美国威斯康辛州森林景观

mba.int <- mba.surf(LIDAR, 300, 300, extend = TRUE)$xyz.est image(mba.int, xaxs = "r", yaxs = "r") # 透视图 persp(mba.int, theta = 90, phi = 20, col = "green3", scale = FALSE, ltheta = -120, shade = 0.75, expand = 10, border = NA, box = FALSE ) z <- mba.int$z
nrz <- nrow(z)
ncz <- ncol(z)
nbcol <- 100
color <- hcl.colors(nbcol)
# Compute the z-value at the facet centres
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# z 分量用颜色表示
facetcol <- cut(zfacet, nbcol)
# png(file="mba.png",res=300,width = 680,height = 680)
par(mar = rep(0, 4))
## 透视图
persp(mba.int,
theta = 90, phi = 20, col = color[facetcol], scale = FALSE,
ltheta = -120, shade = 0.75, expand = 20, border = NA, box = FALSE
)
# 热图
image(mba.int, xaxs = "r", yaxs = "r", col = gray(seq(1, 0, l = 101)))

library(barsurf)
x <- y <- 1:4
f <- function(x, y) x^2 + y^2
z <- outer(x, y, f)
plot3d.bar(, , z)
plot3d.bar(, , volcano)
plot3d.surf(, , volcano)

## 11.6 TikZ 绘图

knitr 提供 Tikz 图形的模版， system.file('misc', 'tikz2pdf.tex', package = 'knitr')tikzDevice 包可以方便的把 R 代码转化为 tikz 代码，然后使用 LaTeX 引擎编译成 PDF 文档，特别地，它很好地支持了图里的数学公式

library(tikzDevice)
tf <- file.path(getwd(), "demo-tikzDevice.tex")

tikz(tf, width = 6, height = 4, pointsize = 30, standAlone = TRUE)
# 绘图的代码，仅支持 Base R Graphics System
source(file = "code/chapter_03/matern.R")
dev.off()

tools::texi2dvi(tf, pdf = T)

system("rm demo-tikzDevice.tex *.log *.aux *.dvi")

system("convert -density 300 -trim demo-tikzDevice.pdf -quality 100 demo-tikzDevice.png")

system("mv demo-tikzDevice.* figures/")
# convert test.svg test.png

# 带有图标题
x <- rnorm(10)
y <- x + rnorm(5, sd = 0.25)
model <- lm(y ~ x)
rsq <- summary(model)$r.squared rsq <- signif(rsq, 4) plot(x, y, main = "Hello \\LaTeX!", xlab = "$x$", ylab = "$y$") abline(model, col = "red") mtext(paste("Linear model:$R^{2}=", rsq, "$"), line = 0.5) legend("bottomright", legend = paste("$y = ", round(coef(model)[2], 3),
"x +", round(coef(model)[1], 3), "$", sep = "" ), bty = "n") plot(x, y, main = "Hello \\LaTeX!", xlab = "$x$", ylab = "$y$") abline(model, col = "red") mtext(paste("Linear model:$R^{2}=", rsq, "$"), line = 0.5) legend("bottomright", legend = paste("$y = ", round(coef(model)[2], 3),
"x +", round(coef(model)[1], 3), "$", sep = "" ), bty = "n") x0 <- 2^(-20:10) nus <- c(0:5, 10, 20) x <- seq(0, 4, length.out = 501) plot(x0, x0^-8, frame.plot = TRUE, # 添加绘图框 log = "xy", # x 和 y 轴都取对数尺度 axes = FALSE, # 去掉坐标轴 xlab = "$u$", ylab = "$\\mathcal{K}_{\\kappa}(u)$", # 设置坐标轴标签 type = "n", # 清除绘图区域的内容 ann = TRUE, # 添加标题 x和y轴标签 panel.first = grid() # 添加背景参考线 ) axis(1, at = 10^seq(from = -8, to = 2, by = 2), labels = paste0("$\\mathsf{10^{", seq(from = -8, to = 2, by = 2), "}}$") ) axis(2, at = 10^seq(from = -8, to = 56, by = 16), labels = paste0("$\\mathsf{10^{", seq(from = -8, to = 56, by = 16), "}}$"), las = 1 ) for (i in seq(length(nus))) { lines(x0, besselK(x0, nu = nus[i]), col = hcl.colors(9)[i], lwd = 2) } legend("topright", legend = paste0("$\\kappa=", rev(nus), "$"), col = hcl.colors(9, rev = T), lwd = 2, cex = 1 ) x <- seq(0, 4, length.out = 501) x <- x[x > 0] plot(x, x, frame.plot = TRUE, ylim = c(1e+0, 1e+20), log = "y", xlab = "$u$", type = "n", yaxt = "n", ylab = "$\\mathcal{K}_{\\kappa}(u)$", ann = TRUE, panel.first = grid() ) axis(2, at = c(1e+0, 1e+05, 1e+10, 1e+15, 1e+20), labels = paste0("$\\mathsf{10^{", seq(from = 0, to = 20, by = 5), "}}$"), las = 1 ) for (i in seq(length(nus))) { lines(x, besselK(x, nu = nus[i]), col = hcl.colors(9)[i], lwd = 2) } legend("topright", legend = paste0("$\\kappa=", rev(nus), "$"), col = hcl.colors(9, rev = T), lwd = 2, cex = 1 ) \usetikzlibrary{arrows} \begin{tikzpicture}[node distance=2cm, auto,>=latex', thick, scale = 0.5] \node (P) {$P$}; \node (B) [right of=P] {$B$}; \node (A) [below of=P] {$A$}; \node (C) [below of=B] {$C$}; \node (P1) [node distance=1.4cm, left of=P, above of=P] {$\hat{P}$}; \draw[->] (P) to node {$f$} (B); \draw[->] (P) to node [swap] {$g$} (A); \draw[->] (A) to node [swap] {$f$} (C); \draw[->] (B) to node {$g$} (C); \draw[->, bend right] (P1) to node [swap] {$\hat{g}$} (A); \draw[->, bend left] (P1) to node {$\hat{f}$} (B); \draw[->, dashed] (P1) to node {$k\$} (P);
\end{tikzpicture}
\begin{tikzpicture}
\begin{scope}[blend group = soft light]
\fill[red!30!white]   ( 90:1.2) circle (2);
\fill[green!30!white] (210:1.2) circle (2);
\fill[blue!30!white]  (330:1.2) circle (2);
\end{scope}
\node at ( 90:2)    {Typography};
\node at ( 210:2)   {Design};
\node at ( 330:2)   {Coding};
\node [font=\Large] {\LaTeX};
\end{tikzpicture}

## 11.7 图形导出

grSoftVersion()
#>                    cairo                   libpng                     jpeg
#>                 "1.16.0"                 "1.6.34"                    "8.0"
#>                  libtiff
#> "LIBTIFF, Version 4.0.9"

capabilities()
#>        jpeg         png        tiff       tcltk         X11        aqua
#>        TRUE        TRUE        TRUE        TRUE       FALSE       FALSE
#>    http/ftp     sockets      libxml        fifo      cledit       iconv
#>        TRUE        TRUE        TRUE        TRUE       FALSE        TRUE
#>         NLS     profmem       cairo         ICU long.double     libcurl
#>        TRUE        TRUE        TRUE        TRUE        TRUE        TRUE

windows cairo_pdf, cairo_ps
pdf svg
postscript png
xfig jpeg
bitmap bmp
pictex tiff

apropos("dev.")
#>  [1] ".Device"             ".Devices"            "dev.capabilities"
#>  [4] "dev.capture"         "dev.control"         "dev.copy"
#>  [7] "dev.copy2eps"        "dev.copy2pdf"        "dev.cur"
#> [10] "dev.flush"           "dev.hold"            "dev.interactive"
#> [13] "dev.list"            "dev.new"             "dev.next"
#> [16] "dev.off"             "dev.prev"            "dev.print"
#> [19] "dev.set"             "dev.size"            "dev2bitmap"
#> [25] "trellis.device"

The Butterfly Affectation: A case study in embedding an external image in an R plot

Improved Importing of Vector Graphics in R

## 11.8 软件信息

sessionInfo()
#> R Under development (unstable) (2019-11-11 r77397)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 8.1 x64 (build 9600)
#>
#> Matrix products: default
#>
#> locale:
#> [1] LC_COLLATE=Chinese (Simplified)_China.936
#> [2] LC_CTYPE=Chinese (Simplified)_China.936
#> [3] LC_MONETARY=Chinese (Simplified)_China.936
#> [4] LC_NUMERIC=C
#> [5] LC_TIME=Chinese (Simplified)_China.936
#>
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base
#>
#> other attached packages:
#> [1] lattice_0.20-38      latex2exp_0.4.0      rsm_2.10
#> [4] MBA_0.0-9            plot3D_1.1.1         scatterplot3d_0.3-41
#> [7] plotrix_3.7-6        barsurf_0.3.1
#>
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.3        highr_0.8         compiler_4.0.0    pdftools_2.3
#>  [5] tools_4.0.0       digest_0.6.22     evaluate_0.14     rlang_0.4.1
#>  [9] filehash_2.4-2    Matrix_1.2-17     magick_2.2        curl_4.2
#> [13] yaml_2.2.0        mvtnorm_1.0-11    xfun_0.11         coda_0.19-3
#> [17] stringr_1.4.0     knitr_1.26        askpass_1.1       grid_4.0.0
#> [21] qpdf_1.1          survival_3.1-7    rmarkdown_1.17    bookdown_0.15
#> [25] multcomp_1.4-10   TH.data_1.0-10    magrittr_1.5      codetools_0.2-16
#> [29] emmeans_1.4.2     htmltools_0.4.0   splines_4.0.0     MASS_7.3-51.4
#> [33] misc3d_0.8-4      tikzDevice_0.12.3 xtable_1.8-4      colorspace_1.4-1
#> [37] tinytex_0.17      sandwich_2.5-1    estimability_1.3  stringi_1.4.3
#> [41] zoo_1.8-6