R语言绘制概率密度图

leengsmile
2016年9月24日

本文源于在高斯混合模型估计中，绘制概率密度图的方法。一般而言，在R语言中，可以使用ggplot2绘制常见的统计图形，本文的方式稍微有所不同，先介绍lattice如何绘制概率密度图，然后再介绍通过ggplot2如何绘制概率密度图。

为了可以重复本文的数据，先设置随机数种子

set.seed(123)

加载magrittr包，以管道形式操作数据流。

require(magrittr)

首先产生如下的一组数据

n <- 1000
mean_s <- c(1, 7)
y <- sample(c("head", "tail"), size = n, replace = TRUE, prob = c(0.25, 0.75))
x <- rnorm(n = 1000, mean = mean_s[1])
tails <- y %in% c("tail")
x[tails] <- rnorm(sum(tails), mean = mean_s[2])

data <- data.frame(x = x, y = y)

这是有服从均值分别为1和7的正态分布随机变量构成的混合变量。在进行数据分析时，可以先查看数据的分布，这里通过lattice包产看数据的概率密度分布

首先引入lattice包

require(lattice)

然后调用densityplot绘制相应的函数

densityplot(~x, data = data, par.settings = list(plot.symbol = list(col = factor(data$y))))

density-plot-1.png

其中par.settings用以将不同混合成分产生的数据用不同颜色的图形表示。par.settings接受由name = value构成的list，可以接收的name可以通过trellis.par.get获取

trellis.par.get() %>% names()

##  [1] "grid.pars"         "fontsize"          "background"       
##  [4] "panel.background"  "clip"              "add.line"         
##  [7] "add.text"          "plot.polygon"      "box.dot"          
## [10] "box.rectangle"     "box.umbrella"      "dot.line"         
## [13] "dot.symbol"        "plot.line"         "plot.symbol"      
## [16] "reference.line"    "strip.background"  "strip.shingle"    
## [19] "strip.border"      "superpose.line"    "superpose.symbol" 
## [22] "superpose.polygon" "regions"           "shade.colors"     
## [25] "axis.line"         "axis.text"         "axis.components"  
## [28] "layout.heights"    "layout.widths"     "box.3d"           
## [31] "par.xlab.text"     "par.ylab.text"     "par.zlab.text"    
## [34] "par.main.text"     "par.sub.text"

所绘制的概率密度图有plot.symbol控制，可以通过一个list参数列表指定。plot.symbol可以接受的参数通过向trellis.par.get传入"plot.symbol"即可。

trellis.par.get("plot.symbol")

## $alpha
## [1] 1
## 
## $cex
## [1] 0.8
## 
## $col
## [1] "#0080ff"
## 
## $font
## [1] 1
## 
## $pch
## [1] 1
## 
## $fill
## [1] "transparent"

在这里通过col指定不同的混合成分指定不同的颜色，在data数据中，不同的混合成分，通过y变量可以获取。因此将data$y赋值给col。

可以通过plot.symbol指定图形的类型

densityplot(~x, data = data, par.settings = list(plot.symbol = list(col = factor(data$y),
                                                                    pch = 25, 
                                                                    alpha = 0.4)
                                                 
                                                 
                                                 )
)

unnamed-chunk-7-1.png

下面在来看看ggplot2中如何绘制类似的图。

首先引入ggplot2

require(ggplot2)

然后通过geom_density来绘制概率密度图

ggplot(data, aes(x = x)) + geom_density(colour = "cadetblue3")

unnamed-chunk-9-1.png

基本的图形已经出来，但是有以下几个问题

1. 密度曲线下方有条直线。
2. 没有原始数据的分布情况。类似于densityplot中的各个数据点，根据这些数据点，可以直观的判断数据的分布情况。

下面依次解决相应的问题

首先去除概率密度图下方的直线，可以使用geom_line，并制定stat = "density"，可以参见[1]

ggplot(data, aes(x = x)) + geom_line(colour = "cadetblue3", stat = "density")

unnamed-chunk-10-1.png

下一步要标记原始数据，在ggplot2中，我觉得更适合用geom_rug来表示原始数据的分布情况。<span style="col:red">2016年9月24日<span>

ggplot(data, aes(x = x)) + geom_line(colour = "cadetblue3", stat = "density") +
    geom_rug(aes(colour = y), sides = "b")

unnamed-chunk-11-1.png

当然，也可以绘制与densityplot类似的图形

ggplot(data, aes(x = x)) + geom_line(colour = "cadetblue3", stat = "density") +
    geom_point(aes(y = 0, colour = y), shape = 25, alpha = 0.4, size = 4) +
    theme(legend.position = "none")

unnamed-chunk-12-1.png

参考文献

[1]. http://stackoverflow.com/questions/21600754/ggplot2-and-geom-density-how-to-remove-baseline

[2]. http://stackoverflow.com/questions/14604435/turning-off-some-legends-in-a-ggplot