复现 sci 顶刊中的 3D 密度函数图

简介

最近在看一些关于贝叶斯深度学习在可靠性方向应用的文章，看到下面这篇文章，发表在可靠性方向顶刊的 ITR 中。

本文主要介绍使用贝叶斯深度学习来预测产品剩余使用寿命，但是传统的深度学习方法只能得到未来时间点的估计值，并不能包含不确定性因素，而使用贝叶斯深度学习，最后可以得到未来时间点预测值的分布。

文中在模拟阶段给出了一个用 matlab 得到的 3D 密度函数图。

主要目的是解释使用该方法预测出未来时间点对应剩余使用寿命（RUL）的分布，对应的点估计和真实值。可以看出，该方法点估计和真实值非常接近，并且还给出对应点的核密度函数。根据这个核密度函数你可以求出 $95\%$​ 的预测区间。

好像内容介绍太多了，本文不是文献解读😁。主要是复现这个图，那正式开始吧！

教程

对于这种 3D 的图形，ggplot2 包以及对应的拓展包好像不是很在行，但是也有一些教程可见：3D Plots with ggplot2 and Plotly。这里我们主要使用 plot3D 包中的 scatter3D 函数进行绘制，当然也可以尝试使用 Scatterplot3d 包。

首先构造一些模拟数据作为例子。这里我们假设每个时间点的密度函数服从正态分布，均值分别为 1：5，标准差都为 1。颜色是自己比较喜欢的几种配色，参考小明的推文：R语言ggplot2画图一套好看的配色以及调整字体的简单小例子。

mean1 = 1:5 # 刻画不同时间对应密度函数的均值
len = 1000
col = c("#02B1e6", "#E81D22", "#F9BC15", "#8015f9", "#20e81d", "#e64602")
x <- seq(-5,10, length = len)
y <- seq(1,5, length = 5)
z = matrix(NA,length(x),5)
for(i in 1:length(mean1)){
  z[,i] = dnorm(x,mean1[i],1)
}

之后使用 scatter3D 函数构建一个空的三维箱子（alpha=0）。这个函数前三个参数分别对应 x，y，z 坐标，bty（boxtype）表示箱子的类型，这里使用类似文献中的形式，之后我们将展示其他形式。phi，theta 主要是对上下角度和左右角度进行调整，ticktype = "detailed"将坐标轴详细展现出来。

如果你对其他参数感兴趣，可以参见这篇教程：Impressive package for 3D and 4D graph - R software and data visualization。

scatter3D(x, rep(1,len), z[,1], bty = "b2",colkey = FALSE,
          phi = 14,theta = 45,
          pch = 18,alpha = 0,ticktype = "detailed", expand =0.5,
          ylim = c(1,5),zlim = c(0,0.5),
          xlab = "RUL(10s)", ylab = "Time(10s)", zlab = "Kernel Distribution"
)

之后就是使用参数内部的add = TRUE 加入各个密度函数图以及点估计。

for(i in 1:length(mean1)){
  scatter3D(x, rep(i,len), z[,i], add = TRUE,type ='l', col = col[i], pch = 16)
  scatter3D(x[which.max(z[,i])], y[i], 0 , add = TRUE, type ='h', col = col[i], pch = 16)
}

这时基本得到了与文献中类似的 3D 密度函数图啦！ 真实值也可以在循环中加入，这里就不再进行啦。

如果你想换一种箱子风格可以调整 bty 参数，例如：bty = "g" , bty = "f"等。

小编有话说

• 这个图其实还少了一个图例，我还没将其进行添加，在此备注下。
• 这个系列应该会继续下去，主要绘制一些在工业工程方向（小编方向）中常用的图。说明下，这个方向绘制的图会比较简单，难点在于模型建立和理论推导，所以小编也不知复现这些图对读者们帮助大不大，主要是自己科研期间的记录。