BSC-贝叶斯-如何从隐含的样本里抽样？

knitr::opts_chunk$set(tidy = TRUE, 
                      warning = FALSE,
                      message = FALSE)
library(showtext) #载入库

lab 9

lab 10

proposal-logistic

观察

rm(list=ls())
# d是概率密度
# p为概率分布函数。pnorm(z)等价于P[x<=z]
# q为pnorm的逆运算，返回给定P下的z
# r生成符合分布的随机数
fun.logis<-function(x){
  return(-dnorm(x,0,1)/dlogis(x,0,1))
}

M.logis<-(-1)*optimise(fun.logis,c(0,2))$objective

#绘制曲线(逻辑分布和正态分布)
curve(dnorm(x,0,1),-5,5,col = "blue",ylim=c(0,0.5),xlab = "x",ylab = "Density",n=200,lwd=2)
curve(M.logis*dlogis(x,0,1),-5,5,add = TRUE,col="red",n=200,lwd=2)
legend("topright",inset = 0.04,c("real","proposal"),pch = c(15:15),col = c("blue","red"))

• 上述值得说明的一点是，由于M未知。所以需要通过$max\frac{f(\theta |y)}{g(\theta )}$求出最大值

• 此外，为了使用optimise内置函数，这里用到了取负数，最大化正数，即最小化负数。

rejection sampling

res.logis<-c()
count.logis<-0
for(i in 1:10000){
  x=rlogis(1)
  if(runif(1)<dnorm(x)/(M.logis*dlogis(x))){
    count.logis=count.logis+1
    res.logis<-c(res.logis,x)
  }
}

real.logis=count.logis/10000
print(paste("logistic proposal实际接受概率为",real.logis))
#实际接受概率

curve(dnorm(x,0,1),-5,5,lwd=2,col ="blue",xlab = "x",ylab = "Density",n=200,main="logit-proposal")
lines(density(res.logis),lwd=2,col="red")
legend("topright",inset = 0.04,c("real","from-proposal"),pch = c(15:15),col = c("blue","red"))

• 10000次的抽样，从$g(\theta )$也就是logistic的分布中抽取样本，然后计算接受概率即dnorm(x)/(M.logis*dlogis(x))

• 最后可视化这10000次抽出的样本的分布密度图如上

• 总计接受了6254个样本，M=1.5957，正好为其倒数，和上课的内容$1/M$一致

proposal-cauchy

观察

rm(list=ls())
# d是概率密度
# p为概率分布函数。pnorm(z)等价于P[x<=z]
# q为pnorm的逆运算，返回给定P下的z
# r生成符合分布的随机数

fun.cauchy<-function(x){
  return(-dnorm(x,0,1)/dcauchy(x,0,1))
}

M.cauchy<-(-1)*optimise(fun.cauchy,c(0,10))$objective

#绘制曲线
curve(dnorm(x,0,1),-5,5,col = "blue",lwd=2,ylim=c(0,0.5),xlab = "x",ylab = "Density",n=200)
curve(M.cauchy*dcauchy(x,0,1),-5,5,lwd=2,add = TRUE,col="red",n=200)
legend("topright",inset = 0.04,c("real","proposal"),pch = c(15:15),col = c("blue","red"))

• cauchy过程同最开始的losigitc，上述值得说明的一点是，由于M未知。所以需要通过$max\frac{f(\theta |y)}{g(\theta )}$求出最大值

• 此外，为了使用optimise内置函数，这里用到了取负数，最大化正数，即最小化负数。

• 绘制的是经M调整过的cauchy分布的密度和norm正态分布的密度

rejection sampling

res.cauchy<-c()
count.cauchy<-0
for(i in 1:10000){
  x=rcauchy(1)
  if(runif(1)<dnorm(x)/(M.cauchy*dcauchy(x))){
    count.cauchy=count.cauchy+1
    res.cauchy<-c(res.cauchy,x)
  }
}

real.cauchy=count.cauchy/10000
print(paste("cauchy proposal实际接受概率为",real.cauchy))
#实际接受概率

curve(dnorm(x,0,1),-6,6,lwd=2,col ="blue",xlab = "x",ylab = "Density",n=200,main="cauchy-proposal")
lines(density(res.cauchy),lwd=2,col="red")
legend("topright",inset = 0.04,c("real","proposal"),pch = c(15:15),col = c("blue","red"))

• 10000次的抽样，从$g(\theta )$也就是cauchy的分布中抽取样本，然后计算接受概率即dnorm(x)/(M.cauchy*dcauchy(x))

• 最后可视化这10000次抽出的样本的分布密度图如上

proposal-t

观察

rm(list=ls())
# d是概率密度
# p为概率分布函数。pnorm(z)等价于P[x<=z]
# q为pnorm的逆运算，返回给定P下的z
# r生成符合分布的随机数

fun.t<-function(x){
  return(-dnorm(x,0,1)/dt(x,5))
}

M.t<-(-1)*optimise(fun.t,c(0,10))$objective

#绘制曲线
curve(dnorm(x,0,1),-5,5,col = "blue",lwd=2,ylim=c(0,0.5),xlab = "x",ylab = "Density",n=200)
curve(M.t*dt(x,5),-5,5,lwd=2,add = TRUE,col="red",n=200)
legend("topright",inset = 0.04,c("real","proposal"),pch = c(15:15),col = c("blue","red"))

• t分布同最开始的losigitc，上述值得说明的一点是，由于M未知。所以需要通过$max\frac{f(\theta |y)}{g(\theta )}$求出最大值

• 此外，为了使用optimise内置函数，这里用到了取负数，最大化正数，即最小化负数。

• 绘制的是经M调整过的t分布的密度和norm正态分布的密度

rejection sampling

res.t<-c()
count.t<-0
for(i in 1:10000){
  x=rt(1,5)
  if(runif(1)<dnorm(x)/(M.t*dt(x,5))){
    count.t=count.t+1
    res.t<-c(res.t,x)
  }
}

real.t=count.t/10000
print(paste("t proposal实际接受概率为",real.t))
#实际接受概率

curve(dnorm(x,0,1),-6,6,lwd=2,col ="blue",xlab = "x",ylab = "Density",n=200,main="t-proposal")
lines(density(res.t),lwd=2,col="red")
legend("topright",inset = 0.04,c("real","proposal"),pch = c(15:15),col = c("blue","red"))

• 10000次的抽样，从$g(\theta )$也就是t的分布中抽取样本，然后计算接受概率即dnorm(x)/(M.t*dt(x,5))

• 最后可视化这10000次抽出的样本的分布密度图如上

proposal-N(1,2)

观察

rm(list=ls())
# d是概率密度
# p为概率分布函数。pnorm(z)等价于P[x<=z]
# q为pnorm的逆运算，返回给定P下的z
# r生成符合分布的随机数

fun.n<-function(x){
  return(-dnorm(x,0,1)/dnorm(x,1,sqrt(2)))
}

M.n<-(-1)*optimise(fun.n,c(0,10))$objective

#绘制曲线
curve(dnorm(x,0,1),-5,5,col = "blue",lwd=2,ylim=c(0,0.5),xlab = "x",ylab = "Density",n=200)
curve(M.n*dnorm(x,1,sqrt(2)),-5,5,lwd=2,add = TRUE,col="red",n=200)
legend("topright",inset = 0.04,c("real","proposal"),pch = c(15:15),col = c("blue","red"))

• N(1,2)分布同最开始的losigitc，上述值得说明的一点是，由于M未知。所以需要通过$max\frac{f(\theta |y)}{g(\theta )}$求出最大值

• 此外，为了使用optimise内置函数，这里用到了取负数，最大化正数，即最小化负数。

• 绘制的是经M调整过的norm(1,2)分布的密度和norm(0,1)分布的密度

rejection sampling

res.n<-c()
count.n<-0
for(i in 1:10000){
  x=rnorm(1,1,sqrt(2))
  if(runif(1)<dnorm(x)/(M.n*dnorm(x,1,sqrt(2)))){
    count.n=count.n+1
    res.n<-c(res.n,x)
  }
}

real.n=count.n/10000
print(paste("N(1,2) proposal实际接受概率为",real.n))
#实际接受概率

curve(dnorm(x,0,1),-6,6,lwd=2,col ="blue",xlab = "x",ylab = "Density",n=200,main="N(1,2)-proposal")
lines(density(res.n),lwd=2,col="red")
legend("topright",inset = 0.04,c("real","proposal"),pch = c(15:15),col = c("blue","red"))

• 10000次的抽样，从$g(\theta )$也就是N(1,2)的分布中抽取样本，然后计算接受概率即dnorm(x)/(M.n*dnorm(x,1,sqrt(2)))

• 最后可视化这10000次抽出的样本的分布密度图如上