【Learning PGM in R】第七章概率混合模型 R代码及纠错

本文链接：https://blog.csdn.net/qq_43749398/article/details/121990938

R语言小白学习语法又学习用法的艰辛历程。

书是《Learning Probabilistic Graphical Model in R》，书里有很适合小白练手的代码展示，但是部分代码有错误（其实还挺多）。本系列会把能改过来的代码记录下来，给一起学习的小伙伴提供帮助~

每一步都有讲解，不想看讲解的小伙伴跳到每个任务的完整代码部分就好啦~

文章目录

代码任务1 混合模型的期望最大化
代码任务2 专家混合
代码任务3 潜在狄利克雷分配LDA
参考文献

代码任务1 混合模型的期望最大化

分步讲解

本段代码首先从两个高斯分布中随机抽样1000和850个样本，其中 $x1\sim N(-3,2)$ ， $x2\sim N(3.1)$ （对应代码段1）。

#生成两个单变元高斯数据集
x1 <-rnorm(1000, -3, 2) 
x2 <-rnorm(850, 3, 1)
#使用函数 hist 画出x1,x2的频率柱状图
hist(c(x1, x2), breaks =100, col ='red')

Alt

图1-1 单变元高斯数据集

然后导入 $m i x t o o l s$ 包（没有这个包的话，第一行是下载），运用 $n o r m a l m i x E M$ 函数进行结构学习，其中 $x 1$ ， $x 2$ 的比例各占50%（ $l a m b d a = 0.5$ ），并且认为采样结果是由两种模型混合而成（ $k = 2$ ），实际的混合数量也是这样。

#install.packages("mixtools")
library(mixtools)
model <-normalmixEM(c(x1, x2), lambda =0.5, k =2)
#View(model)

观察model的结果（ $v i e w (m o d e l)$ ），我们得到了以下参数。
Alt

图1-2 model的参数

可以看到，混合比例是54.1%和45.9%。这个和书上不一样很正常，因为 $x 1$ ， $x 2$ 就是随机生成的，但是都和初始设定的 $l a m b d a$ 相似。 $\mu$ 的结果是-3.00和3.02，可以说是非常接近初始值了，这是一个很完美的估计结果。我们可以画出采样柱状图和EM算法估计出的高斯分布曲线来直接体现估计的准确性。

#画出柱状图并叠加高斯分布线
hist(c(x1, x2), breaks =100, col ="red", freq = F, ylim =c(0, 0.4))
x = seq(-10,10,0.01)#错误点1
lines(x, dnorm(x, model$mu[1], model$sigma[1]), col ="blue")
lines(x, dnorm(x, model$mu[2], model$sigma[2]), col ="green")

错误1.1：没有定义变量 $x$

书中未定义 $x$ ，所以后续两行 $l i n e s ()$ 函数无法执行。观察代码发现是要画高斯分布曲线，那么只需要取平均分布的 $x$ 坐标即可。最终图片成功显示。
Alt

图1-3 采样柱状图与估计的高斯分布曲线
当然，这是没有添加比例的情况，我们添加比例重新作图。

#添加比例
hist(c(x1, x2), breaks =100, col ="red", freq = F)

lines(x, 
      model$lambda[1] *dnorm(x, model$mu[1], model$sigma[1]) +model$lambda[2]*
        dnorm(x, model$mu[2], model$sigma[2]), 
      lwd =3)

Alt

图1-4 采样柱状图与添加比例的估计的高斯分布曲线

还记得我们在进行参数估计时设置了簇的数量 $k = 2$ ，因为这和实际情况相同，所以最后的结果不错，丹如果模型的数量预设错误，结果会不太理想。所以粗的数量很重要，一般我们会多设几次不同的 $k$ ，寻找最理想的结果。下面运行几个簇与实际情况不相符的代码。

model <-normalmixEM(c(x1,x2), lambda=0.5, k=3) 
model <-normalmixEM(c(x1,x2), lambda=0.5, k=4)

Alt

图1-5 k=2/3情况下的估计结果（均无法收敛）
第一次的结果，两种情况都不能收敛。说明$k$如果设置错误，是有无法收敛（not convergent）的可能性的。

Alt

图1-6 k=3/4情况下的估计结果（k=3收敛）
最后一次，k=3终于收敛了，我们来进行下一步。

补充1.2：画出k=3的图

#画出k=3的图
hist(c(x1, x2), breaks =200, col ="red", freq = F)

lines(x, model1$lambda[1]*dnorm(x, model1$mu[1], model1$sigma[1]) 
        +model1$lambda[2]*dnorm(x, model1$mu[2], model1$sigma[2])
        +model1$lambda[3]*dnorm(x, model1$mu[3], model1$sigma[3]), 
      lwd =3,col="yellow")

结果如图。

Alt

图1-7 k=3情况下的估计结果图像（k=3收敛）
最后结果好像不是很明显，观察参数发现，第三个分布的占比非常之小（0.16%），几乎可以忽略不计，所以$k=3$和$k=2$的差距并不明显。书中给出了当$k$大于实际情况时EM算法的做法猜测：期望最大化算法好像是把较大的簇拆分成了多个高斯分布，进而得到多于实际情况的簇。

完整代码

#生成两个单变元高斯数据集
x1 <-rnorm(1000, -3, 2) 
x2 <-rnorm(850, 3, 1)
#使用函数 hist 画出结果
hist(c(x1, x2), breaks =100, col ='red')
#install.packages("mixtools")
library(mixtools)
model <-normalmixEM(c(x1, x2), lambda =0.5, k =2)
#View(model)
#画出柱状图并叠加高斯分布线
hist(c(x1, x2), breaks =100, col ="red", freq = F, ylim =c(0, 0.4))
x = seq(-10,10,0.01)
lines(x, dnorm(x, model$mu[1], model$sigma[1]), col ="blue")
lines(x, dnorm(x, model$mu[2], model$sigma[2]), col ="green")
#添加比例
hist(c(x1, x2), breaks =100, col ="red", freq = F)

lines(x, 
      model$lambda[1] *dnorm(x, model$mu[1], model$sigma[1]) +model$lambda[2]*
        dnorm(x, model$mu[2], model$sigma[2]), 
      lwd =3)
#改变簇的数量
model1 <-normalmixEM(c(x1,x2), lambda=0.5, k=3) 
model2 <-normalmixEM(c(x1,x2), lambda=0.5, k=4) 
#画出k=3的图
hist(c(x1, x2), breaks =200, col ="red", freq = F)

lines(x, model1$lambda[1]*dnorm(x, model1$mu[1], model1$sigma[1]) 
        +model1$lambda[2]*dnorm(x, model1$mu[2], model1$sigma[2])
        +model1$lambda[3]*dnorm(x, model1$mu[3], model1$sigma[3]), 
      lwd =3, col ="yellow")

代码任务2 专家混合

分步讲解

#生成示例数据集
x1 =runif(40, 0, 10) 
x2 =runif(40, 10, 20) 
e1 =rnorm(20, 0, 2) 
e2 =rnorm(20, 0, 3) 
y1 =1 +2.5 *x1 +e1 
y2 =35 +- 1.5 *x2 +e2 
#作图
xx =c(x1, x2) 
yy =c(y1, y2)

代码一共生成了两个服从平均分布的随机变量数据集合 $x 1, x 2$ 和两个服从正态分布的随机变量数据集合 $e 1, e 2$ 。 $y 1, y 2$ 是两种分布的混合映射（mapping）。画出结果，进行简单的线性回归。
Alt

图2-1 线性回归结果（不理想）
### 补充2.1：画出线性回归结果

plot(xx,yy)
m =lm(yy ~xx)
xx =seq(0, 40, 1) 
lines(xx, xx *m$coefficients[2]+m$coefficients[1], col =2, lw =2)

可以看到线性回归明显不能刻画数据集的行为。它几乎不能捕捉数据趋势，只是数据集合大致平均。

错误2.2：库包名称拼写错误

库包名称是 $h m e E M$
在这里插入图片描述

图2-2 库包名称拼写错误

补充2.3：进行专家混合

书上这一部分的代码没有专家混合的部分(╯▔皿▔)╯，所以笔者参照 $h e l p (h m e E M)$ 写出了下段代码。
首先这是例子代码。

#Examples源代码
## EM output for NOdata. 
data(NOdata)
attach(NOdata)
set.seed(100)
em.out <- regmixEM(Equivalence, NO)
hme.out <- hmeEM(Equivalence, NO, beta = em.out$beta)
hme.out[3:7]

我们照着上面的模仿。

d <- list(x,y)
d1 <- as.data.frame(d)#创建数据框
colnames(d1) <- c("x","y")#修改数据框列名
attach(d1)
#set.seed(100)
em.out <- regmixEM(x, y)
hme.out <- hmeEM(x, y, beta = em.out$beta,k=2)
hme.out[3:7]

$N O d a t$ 的类型是数据框，所以我们把我们的 $x, y$ 放在数据框变量 $d 1$ 中，并修改数据框的列名。另外，笔者在写这段代码时，把 $x, y$ 的位置调换后， $h m e E M$ 一直无法收敛。运行成功后，我们得到了预期的专家混合模型数据。

完整代码

#生成示例数据集
x1 =runif(40, 0, 10) 
x2 =runif(40, 10, 20) 
e1 =rnorm(20, 0, 2) 
e2 =rnorm(20, 0, 3) 
y1 =1 +2.5 *x1 +e1 
y2 =35 +- 1.5 *x2 +e2 
x =c(x1, x2) 
y =c(y1, y2)
plot(x,y)
m =lm(y ~x)
xx =seq(0, 40, 1) 
lines(xx, xx *m$coefficients[2]+m$coefficients[1], col =2, lw =2)
library(mixtools)
d <- list(x,y)
d1 <- as.data.frame(d)
colnames(d1) <- c("x","y")
attach(d1)
#set.seed(100)
em.out <- regmixEM(x, y)
hme.out <- hmeEM(x,y, beta = em.out$beta,k=2)
hme.out[3:7]

最后书上的图片博主不知道怎么画（不知道是哪个函数），如果有突破，会及时更新。
Alt

图2-3 专家混合模型估计结果

代码任务3 潜在狄利克雷分配LDA

分步讲解

错误3.1：库包名称拼写错误

这段代码需要用到 $t o p i c m o d e l s$ 和 $R t e x t T o o l s$ 两个包。但是 $R t e x t T o o l s$ 已经修改为 $R T e x t T o o l s$ ，需要安装的小伙伴，可以把第一二行代码注释掉。安装并导入库包后，纽约时报的数据就可以正常导入了。

install.packages("topicmodels")
install.packages("RTextTools")
library(RTextTools)
library(topicmodels)
data(NYTimes)

错误3.2：随机采样语法错误

书上的采样语法会显示没有 $s a m p l e s ()$ 函数。笔者找到了替代的语句，成功采样。

data <-NYTimes[sample(nrow(NYTimes), 1000),]

错误3.3：language属性值错误

将 $" e n g l i s h, "$ 里的逗号去掉。

matrix <-create_matrix(cbind(as.vector(data$Title),as.vector(data$Subject)), language="english",removeNumbers=TRUE, stemWords=TRUE)
k <-length(unique(data$Topic.Code))
lda <-LDA(matrix, k)
plot(colSums(lda@gamma)/nrow(lda@gamma), t ="h")

$m a t r i x$ 获取了文本标题和文本的话题， $k$ 是唯一主题代码。画出拟合结果中每个主题占的比例。
Alt

图3-1 LDA 拟合结果—各主题占比

补充3.4：随着带有两个以上的主题概率的提高文档数量的变化

文中只给出了图，没有代码。

prob=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)
li= c()
for (p in prob){
  temp=c(sum(sapply(1:nrow(lda@gamma), function(i) sum(lda@gamma[i, ] >p) >1)))
  li <-append(li,temp)
}
plot(prob,li,col=1,type="b",xlab="带有两个以上的概率大于x的主题",ylab="文档数量")

Alt

图3-2 随着带有两个以上的主题概率的提高文档数量的变化

完整代码

#安装并导入库包
install.packages("topicmodels")
install.packages("RTextTools")
library(RTextTools)
library(topicmodels)
#载入《纽约时报数据库》
data(NYTimes)
#采集1000个文本
data <-NYTimes[sample(nrow(NYTimes), 1000),]
#创建矩阵，cbind () 把矩阵横向合并成一个大矩阵（列方式），包含标题和主题
matrix <-create_matrix(cbind(as.vector(data$Title),as.vector(data$Subject)), language="english",removeNumbers=TRUE, stemWords=TRUE)
#使用unique函数得到出现过的主题唯一编码序列。
k <-length(unique(data$Topic.Code))
#运行带有变分期望最大化的学习算法
lda <-LDA(matrix, k)
plot(colSums(lda@gamma), t ="h")
plot(colSums(lda@gamma)/nrow(lda@gamma), t ="h")

sum(sapply(1:nrow(lda@gamma), function(i) sum(lda@gamma[i, ] >0.1) >1))
prob=c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)
li= c()
for (p in prob){
  temp=c(sum(sapply(1:nrow(lda@gamma), function(i) sum(lda@gamma[i, ] >p) >1)))
  li <-append(li,temp)
}
plot(prob,li,col=1,type="b",xlab="带有两个以上的概率大于x的主题",ylab="文档数量")