概率图模型推断算法简单原理+实现

最新推荐文章于 2022-09-27 18:27:04 发布

桃花折江山

最新推荐文章于 2022-09-27 18:27:04 发布

阅读量911

点赞数 1

文章标签：算法图论深度学习概率论

本文链接：https://blog.csdn.net/weixin_51290449/article/details/117330578

版权

本文是一个记录薄，记录一下我关于概率图模型的推断算法的理解

利用连接树算法进行精确推断

联合树算法

library("gRbase")
library("gRain")
library("Rgraphviz")
dag <- dag(~ F+ E:F + C:F + A:C + D:E+B:A:D)
plot(dag)
daG1 <- dag(~ a:b:c + b:c + c:d)
#dag()表明生成有向图，冒号后面的是父节点，冒号的最前面的是子节点
plot(daG1)

下面给出条件概率分布表，条件概率分布表也就对应着图结构，我们给出如dag那样的图结构。把条件概率分布表用compileCPT()函数进行编译

val=c("true","false")
F = cptable(~F, values=c(10,90),levels=val)
C = cptable(~C|F, values=c(10,90,20,80),levels=val)#矩阵默认按列填充
E = cptable(~E|F, values=c(50,50,30,70),levels=val)
A = cptable(~A|C, values=c(50,50,70,30),levels=val)
D = cptable(~D|E, values=c(60,40,70,30),levels=val)
B = cptable(~B|A:D, values=c(60,40,70,30,20,80,10,90),levels=val)
F
C
B
plist = compileCPT(list(F,E,C,A,D,B))#编译条件概率表
plist
print(plist$F)
print(plist$B)

继续利用grain()函数把图结构转换成连接树算法需要的图结构三角化的原理是：指定一个节点的顺序，然后查看与这个节点相连的两个节点是不是属于一个三角形区域，如果不属于的话，那么他们之间是否存在边，如果不存在那么增加边，添加边即可。具体来说，按照节点顺序为FECADB这样的节点顺序，F相连的两个节点为EC，首先看EC节点是不是属于一个三角区域，不是，且EC之间没有边，那么就在EC增加边。

jtree=grain(plist)#调用连接树算法
plot(jtree)#这个就是三角化的图
jtree

把三角化的图转换成团树，在团树上进行推断

#计算F节点的边际分布
querygrain(jtree, nodes=c("F"), type="marginal")
#计算B节点的边际分布
querygrain(jtree, nodes=c("B"), type="marginal")
#计算C的边际分布
querygrain(jtree, nodes=c("C"), type="marginal")
#对A,B的联合概率分布进行推断，联合概率分布中行与列是对等的关系
querygrain(jtree, nodes=c("A","B"), type="joint")
#对ABC的联合概率分布进行计算
querygrain(jtree, nodes=c("A","B","C"), type="joint")

上面我们主要是求边际分布，边际分布很重要，边际分布是求条件分布的关键，下面我们利用连接树算法求条件分布

#首先设置证据节点F=TRUE，在F已知的前提下求目标变量的边际分布，这个时候的边际分布也叫做后验概率分布，这是我们推断中最关心的一个问题
jtree2 = setEvidence(jtree, evidence=list(F="true"))
#对比原来的边际分布也叫做先验分布和现在在F已知的条件下的后验分布
querygrain(jtree, nodes=c("F"), type="marginal")
querygrain(jtree2, nodes=c("F"), type="marginal")

querygrain(jtree, nodes=c("A"), type="marginal")
querygrain(jtree2, nodes=c("A"), type="marginal")

querygrain(jtree, nodes=c("B"), type="marginal")
querygrain(jtree2, nodes=c("B"), type="marginal")

当我们给更多的证据时，会对目标节点的后验分布又会变的不一样

jtree3 = setEvidence(jtree, evidence=list(F="true",A="false"))

querygrain(jtree, nodes=c("C"), type="marginal")
querygrain(jtree2, nodes=c("C"), type="marginal")
querygrain(jtree3, nodes=c("C"), type="marginal")

下面看一下包gRbase里的数据，分别给它有向图和无向图的结构，利用querygrain()函数进行推断

data(lizard, package="gRbase")
daG <- dag(~species + height:species + diam:species)
plot(daG)
uG  <- ug(~height:species + diam:species)
plot(uG)
bn.uG   <- grain(uG, data=lizard)
querygrain(bn.uG,nodes="height",type="marginal")
bn.daG  <- grain(daG, data=lizard)
#无向图也可以进行边际分布的计算
querygrain(bn.daG,nodes="height",type="marginal")

近似推断算法

当我们进行边际分布的推断时，分母部分是证据变量的联合分布，当证据变量只有一个时，分母就是证据变量的边际分布，不管是求证据变量的联合分布还是边际分布都需要进行要么对无关变量求和，要么对无关变量进行积分的操作，在连续变量下，经常会遇到积分不容易求，计算量非常大的情况，那么在这个时候就引入了我们的近似推断算法。因为我们有时计算某个目标变量的分布，并不是单纯的计算它的分布，我们往往是要利用这个分布，计算某个函数在这个分布下的期望，那么我们可以利用近似的思想，先从目标分布中抽取样本，把样本带入函数中求函数的平均值，也叫做样本平均值，利用大树定理可知当样本容量增大时，样本均值以概率1收敛于数学期望。所以问题就演变成怎么从目标分布中抽取样本的问题（注意这里的目标分布可以是未归一化的分布）

直接采样法

直接采样法适用对象：离散分布（比较简单一般不讨论），连续分布密度函数已知且分布函数可求，且分布函数可以求出反函数。

从离散分布中随机采样

常见的离散分布有二项分布，泊松分布，几何分布最常用的应该是二项分布，多项分布

x=rbinom(100,10,0.3)
hist(x,prob=T,main=paste("n =",10))
#下面展示怎么从多项分布中抽样，可以直接调用sample()函数
x1=sample(1:6,100,replace=TRUE)
hist(x1,prob=T,main=paste("n =",100))
x1=sample(1:6,4,replace=TRUE,prob=c(0.1,0.3,0.2,0.2,0.1,0.3))
x1

从连续分布中利用均匀分布抽样

利用0-1均匀分布从指数分布中抽样从下面的我们采样得到数据的直方图和同样参数的指数分布图可以明显的看出利用来自均匀分布的随机数的变换抽取的随机数是服从我们的目标分布指数分布的

x2=runif(20000)
inv_h=function(x,lambda) {-(1/lambda)*log(1-x2)}#这个inv_h服从的就是指数分布
#画出直方图，且纵轴表示的是概率
hist(inv_h(x2,2),freq=FALSE,breaks=100)
t=seq(0,4,0.01)
lines(t,dexp(t,2),lw=2)

间接采样法

拒绝采样法

原理：假设我们要从p(x)分布中取样，直接取样不好取的时候我们引入一个可以直接抽样的分布，称为建议分布q(X),并且有一个数k,一定满足kq(x)>p(x). 按照q(x)进行抽样，假设得到的结果是x’,再按照p(x’)/kq(x’)的比例随机决定是否接受抽取的样本x’

#首先定义建议分布和目标分布
q=function(x) {dnorm(x,0,0.5)}  #定义建议分布
 
rq=function(x) {rnorm(1,0,0.5)}  #从建议分布中随机取出一个数

p=function(x) {0.6*dnorm(x,0,0.1)+0.4*dnorm(x,0.6,0.2)} #定义目标分布 
#然后写出拒接采样的算法
rejection=function(N,k,p,q,rq) {
  #创建两个向量用来接受收取的样本和是否接受这个样本 
  accept=logical(N)
  x=numeric(N)

  for(i in 1:N) {
    z0=rq()
    u0=runif(1,0,1)
    if (u0<p(z0)/k*q(z0)) {
      accept[i]=TRUE
    }
    else {
      accept[i]=FALSE
    }
    x[i]=z0
  }
  return(data.frame(x=x,accept=accept))
  }

set.seed(600)
k=3.1
x3=rejection(100,k,p,q,rq)#接受返回值的
sum(x3$accept)
#增加采样数
x4=rejection(5000,k,p,q,rq)
sum(x4$accept)
x5=rejection(50000,k,p,q,rq)
sum(x5$accept)
t1=seq(-2,2,0.01)
#展示不同抽样次数，产生的样本与目标分布之间的关系
hist(x3$x[x3$accept],freq=F,breaks=200,col="grey")
lines(t1,p(t1),col=2,lwd=2)

hist(x4$x[x4$accept],freq=F,breaks=200,col="grey")
lines(t1,p(t1),col=2,lwd=2)

hist(x5$x[x5$accept],freq=F,breaks=200,col="grey")
lines(t1,p(t1),col=2,lwd=2)

所以当k较大时，需要采样很多样本才能使保留下来的样本近似服从目标分布，所以拒绝采样效率可能不高，而且找到一个和目标分布比较接近的提议分布也比较困难。

重要性采样

当我们的目标就是计算目标分布p(x)下函数f(x)的期望，实际上抽取的样本并不需要严格服从目标分布p(x),也可以通过从建议分布中采样并估计期望，我们此时不要求建议分布与目标分布比较接近。

#定义重要性采样算法
#参数说明
#f:是想要求的期望的函数，N是抽取的样本数目，P是目标分布，q是建议分布
#我们从下面的例子中计算平均值，
#为了方便看出计算的效果，使用的函数为f(x)=x

importance=function(N,f,p,q,rq) {
  x=sapply(1:N,rq)#从提议分布中采样
  A=sum((f(x)*p(x)/q(x)))#分子
  B=sum(p(x)/q(x))#分母
  return(A/B)
}
#例子1：想要从混合高斯分布下计算函数的期望，引入建议分布高斯分布
set.seed(600)
p1=function(x) {0.6*dnorm(x,0,0.1)+0.4*dnorm(x,0.6,0.2)}
q1=function(x) {dnorm(x,0,0.5)}
rq1=function(x) {rnorm(1,0,0.5)}
f=function(x) {x}
plot(t1,p1(t1),type="l",col="red")
lines(t1,q1(t1),type="l",col="black")
print(importance(1000,f,p1,q1,rq1))
print(importance(10000,f,p1,q1,rq1))
print(importance(50000,f,p1,q1,rq1))

#下面我们看一下样本均值怎么收敛到真实均值的
t2=seq(1000,50000,500)
#1、先看实验1
x6=sapply(t2,function(i) {importance(i,f,p1,q1,rq1)})
plot(x6,type="l")
abline(h=0.24,col="red")#添加水平线

以上我在拒绝采样和重要性采样中举的例子都是一维的，当在高维空间中拒绝采样和重要性采样的效率随空间维数的增加而指数降低，此时马尔科夫蒙塔卡洛方法（MCMC）是一种更好的采样方法，可以很容易对高维变量进行采样。

MCMC采样

MCMC包括MH算法和吉布斯采样（Gibbs Sample）核心思想是把采样过程看做一个马尔科夫链，也就是说t+1时刻抽取的样本与t时刻抽取的样本有关，拿天气状况举个例子：今天的天气为晴天，明天的天气状况只依赖于今天的天气状况和状况之间的转移分布，如果是离散变量这个状态之间的转移分布是状态转移矩阵，如果是连续变量，那么状态转移分布为转移核函数。一个不可约非周期其正常返的马尔科夫链一定有平稳分布，如果这个平稳分布是我们的目标分布p(x),那么在状态平稳时抽取的样本就是服从p(x)的，所以目标是从p(x)中抽取样本，那怎么抽？在拒绝采样的时候样本其实是从建议分布中抽的，只是采取策略让接受的样本近似从p(x)中抽取，但是MCMC方法就不一样了，MCMC最后抽取的样本就是从目标分布p(x)中抽取的，我们只需要构造平稳分布为p(x)的马尔科夫链即可，也就是找到合适的状态转移矩阵或者是转移核函数q(x’|x)。

MH算法

#例子：假设目标分布为N(3,5),建议函数采用条件概率分布p(x'|x),定义为多元正态分布，均值为x,协方差矩阵为常数矩阵
#即新的状态 x'~N(x,R),R为标准差，为常数即可，这里给的是2
N=10000
x7=c()
x0=100
x7[1]=x0
for (i in 2:N) {
  u=runif(1)
  y=rnorm(1,x7[i-1],2)
  A=min(1,dnorm(y,3,5)/dnorm(x7[i-1],3,5))
  if(u<A) {
    x7[i]=y
  }
  else {
    x7[i]=x7[i-1]
  }
}

hist(x7[700:N],freq=F)
curve(dnorm(x,3,5),add=TRUE)

acf(x7[700:N])

z=x7[seq(700,N,by=40)]
acf(z)
hist(z,freq=F)
curve(dnorm(x,3,5),add=TRUE)

ks.test(z,"pnorm",3,5)#检验选取样本z确实是来自我们的正态分布的样本

2、metropolis选择的另一个特例是随机游走metropolis算法，该算法是假设t时刻的样本为x,t+1时刻的抽取的样本是在原x的基础上加上一个随机游走项，这个随机游走可以是多种形式的随机游走，最简单那的是在一个均匀分布范围内进行随机游走，或者是以标准正态变量的比例进行随机游走，下面对这两种形式的随机游走metropolis算法进行实现

#例子
#从贝塔分布beta(2,7)中抽取样本，初始点为0.1,随机游走在u(-0.2,0.2)
N1 =30000 #抽样个数
x8 =c()
x01 =0.1 # 初始值
x8[1]=x01
k = 0 #k表示拒绝转移的次数
for (i in 2:N) {
  y = x8[i-1]+runif(1,-0.2,0.2)
  u=runif(1)
  if(0<y&y<1){
    A=min(1,dbeta(y,2,7)/dbeta(x8[i-1],2,7))
    if (u<A)
      x8[i]= y
    else {
      x8[i]=x8[i-1]
      k =k + 1
    }
  }       
  else
    x8[i]=x01
}

plot(x8,type="l")
acf(x8)
z=x8[seq(5000,N,by=20)]
ks.test(z,"pbeta",2,7)
ks.test(x8[1000:N],"pbeta",2,7)

#例子：随机样本x1,x2,,,,xn来自均值为u,方差为sigma^2的正态分布，现在已知观测数据
#为x1,x2,...xn，u的先验分布为u(-10,10),sigma~exp(1)
#目的：估计u和sigma

#下面利用metrop()函数来进行估计
#下面定义的f是后验分布的对数似然函数
library("mcmc")
f=function(theta,x) {
  mu=theta[1]
  sigma=theta[2]
  n=length(x)
  if(abs(mu)<=10&&sigma>0) {
    y=-n*log(sigma)-sigma-sum((x-mu)^2)/(2*sigma^2)
  }
  else {
    y=-Inf
  }
  return(y)
}

n2=100
x9=rnorm(n2,5,2)
N2=1000
theta.int=c(mean(x9),sd(x9))
out=metrop(f,theta.int,nbatch=N2,scale=0.2,x=x9)
theta=out$batch
theta
acf(theta[10:N2,1])
mean(theta[10:N2,1])
mean(theta[10:N2,2])
mean(theta[seq(1,N2,by=20)])

Gibbs采样算法

Gibbs抽样用于对多元变量联合分布的抽样和估计，基本做法是：从联合概率分布定义满条件概率分布，以此对满条件概率分布进行抽样，得到样本序列，这样的抽样过程是在一个马尔科夫链上的随机游走，每一个样本对应着马尔科夫链的状态，平稳分布就是目标分布，在Gibbs抽样中我们的目标分布往往是联合分布，燃烧期之后的样本就是联合分布的随机样本。

Gibbs抽样是MCMC的一种，故Gibbs抽样仍然是利用建议分布，只不过这里的建议分布为目标分布的满条件概率分布下面看一个具体的例子用吉布斯抽样从以下二维正态分布中随机抽取样本

#统计学习方法p377页例子的实现
library("MASS")
N2=5000
burn.in=2500
X=matrix(0,N2,2)
mu1=0
mu2=0
sigma1=1
sigma2=1
rho=0.5#相关系数
s1.c=sqrt(1-rho^2)*sigma1
s2.c=sqrt(1-rho^2)*sigma2
X[1,]=c(mu1,mu2)#初始化
for (i in 2:N2) {
  #首先在x2的条件下抽x1
  x2=X[i-1,2]
  m1.c=mu1+rho*(x2-mu2)*sigma1/sigma2
  X[i,1]=rnorm(1,m1.c,s1.c)
  x1=X[i,1]
  m2.c=mu2+rho*(x1-mu1)*sigma2/sigma1
  X[i,2]=rnorm(1,m2.c,s2.c)
}
b=burn.in+1
x.mcmc=X[b:N2,]
#下面画出二维正态分布的 密度曲线，画出等高线
bigauss.estimate=kde2d(x.mcmc[,1],x.mcmc[,2],n=1000)

contour(bigauss.estimate,nlevels=15,lty=2)

桃花折江山

关注

1
点赞
踩
1

收藏

觉得还不错? 一键收藏
0
评论
概率图模型推断算法简单原理+实现

本文是一个记录薄，记录一下我关于概率图模型的推断算法的理解利用连接树算法进行精确推断联合树算法library("gRbase")library("gRain")library("Rgraphviz")dag <- dag(~ F+ E:F + C:F + A:C + D:E+B:A:D)plot(dag)daG1 <- dag(~ a:b:c + b:c + c:d)#dag()表明生成有向图，冒号后面的是父节点，冒号的最前面的是子节点plot(daG1)下面给出条件概率
复制链接

扫一扫