西部数码网站管理助手 xp,亳州市城乡建设局网站,域名是什么结构称为域名空间,不会编程可以做网站吗全文链接#xff1a;http://tecdat.cn/?p23785 最近#xff0c;我们使用贝叶斯非参数#xff08;BNP#xff09;混合模型进行马尔科夫链蒙特卡洛#xff08;MCMC#xff09;推断#xff08;点击文末“阅读原文”获取完整代码数据#xff09;。 概述 相关视频 在这篇文… 全文链接http://tecdat.cn/?p23785 最近我们使用贝叶斯非参数BNP混合模型进行马尔科夫链蒙特卡洛MCMC推断点击文末“阅读原文”获取完整代码数据。 概述 相关视频 在这篇文章中我们通过展示如何使用具有不同内核的非参数混合模型进行密度估计。在后面的文章中我们将采用参数化的广义线性混合模型并展示如何切换到非参数化的随机效应表示避免了正态分布的随机效应假设。 使用Dirichlet Process Mixture模型进行基本密度估计 提供了通过Dirichlet过程混合DPM模型进行非参数密度估计的机制Ferguson, 1974; Lo, 1984; Escobar, 1994; Escobar and West, 1995。对于一个独立和相同分布的样本 该模型的形式为 这个模型实现是灵活的运行任意核的混合。, 可以是共轭的也可以是不共轭的也是任意的基度量 . 在共轭核/基数测量对的情况下能够检测共轭的存在并利用它来提高采样器的性能。 为了说明这些能力我们考虑对R中提供的Faithful火山数据集的喷发间隔时间的概率密度函数进行估计。 data(faithful) 观测值 对应于数据框架的第二列而 . 使用CRP表示法拟合高斯_location-scale_ 分布混合分布 模型说明 我们首先考虑用混合正态分布的_location-scale_Dirichlet过程s来拟合转换后的数据 其中 对应的是正态-逆伽马分布。这个模型可以解释为提供一个贝叶斯版本的核密度估计 用于使用高斯核和自适应带宽。在数据的原始尺度上这可以转化为一个自适应的对数高斯核密度估计。 引入辅助变量表明混合的哪个成分产生了每个观测值并对随机量进行积分我们得到模型的CRP表示Blackwell and MacQueen, 1973。 其中 是向量中唯一值的数量是第个唯一值在中出现的次数。这个说明清楚地表明每个观测值都属于最多正态分布聚类中的任何一个并且CRP分布与分区结构的先验分布相对应。 这个模型的说明是这样的 y\[i\] ~ dnorm(mu\[i\], var s2\[i\])mu\[i\] - muTilde\[xi\[i\]\]s2\[i\] - s2Tilde\[xi\[i\]\]xi\[1:n\] ~ dCRP(alpha, size n)muTilde\[i\] ~ dnorm(0, var s2Tilde\[i\])s2Tilde\[i\] ~ dinvgamma(2, 1)alpha ~ dgamma(1, 1) 请注意在模型代码中参数向量muTilde和s2Tilde的长度被设置为.我们这样做是因为目前的实现要求提前设置参数向量的长度并且不允许它们的数量在迭代之间变化。因此如果我们要确保算法总是按预期执行我们需要在最坏的情况下工作即有多少个成分就有多少个观测值的情况。但它的效率也有点低无论是在内存需求方面当 规模大时需要维护大量未占用的成分还是在计算负担方面每次迭代都需要更新大量不需要后验推理的参数。当我们在下面使用伽马分布的混合时我们将展示一个能提高效率的计算捷径。 还需要注意的是的值控制着我们先验预期的成分数量的值越大对应于数据占据的成分数量越多。因此通过指定一个先验值我们为模型增加了灵活性。对Gamma先验的特殊选择允许使用数据增强方案从相应的全条件分布中有效取样。也可以选择其他的先验在这种情况下这个参数的默认采样是一个自适应的随机游走Metropolis-Hastings算法。 运行MCMC算法 下面的代码设置了数据和常数初始化了参数定义了模型对象并建立和运行了MCMC算法。默认采样器是一个折叠的吉布斯采样器Neal, 2000。 # 模型数据
y standlFaithful
# 模型常量
n length(standlFaithful))
# 参数初始化
list(xi sample(1:10, sizen, replaceTRUE),
# 创建和编译模型
Model(code, data, inits, consts)
##定义模型...
##建立模型...
##设置数据和初始值...
##在模型上运行计算随后的任何错误报告可能只是反映了模型变量的缺失值...
##检查模型的大小和尺寸......。
##模型构建完成。
## 编译完成。
#MCMC的配置、创建和编译
MCMC(conf)
## 编译......这可能需要一分钟
## 编译完成。 我们可以从参数的后验分布中提取样本并创建痕迹图、直方图和任何其他感兴趣的总结。例如对于参数我们有。 # 参数的痕迹图
ts.plot(samples\[ , alpha\], xlab iteration, ylab expression(alpha)) # 后验直方图
hist(samples\[ , alpha\], xlab expression(alpha), main , ylab Frequency) 在这个模型下对于一个新的观察后验预测分布是最佳密度估计在平方误差损失下。这个估计的样本可以很容易地从我们的MCMC产生的样本中计算出来。 # 参数的后验样本samples\[, alpha\]
# 平均值的后验样本samples\[, grep(muTilde, colnames(samples))\] # 聚类平均数的后验样本。
# 集群方差的后验样本
samples\[, grep(s2Tilde, colnames(samples))\] # 聚类成员的后验样本。
# 聚类成员关系的后验样本
samples \[, grep(xi, colnames(samples))\] # 聚类成员的后验样本。hist(y, freq FALSE,xlab 标准化对数尺度上的等待时间)
##对标准化对数网格的密度进行点式估计 然而回顾一下这是对等待时间的对数的密度估计。为了获得原始尺度上的密度我们需要对内核进行适当的转换。 standlGrid*sd(lFaithful) mean(lFaithful) # 对数尺度上的网格hist(faithful$waiting, freq FALSE 无论是哪种情况都有明显的证据表明数据中的等待时间有两个组成部分。 生成混合分布的样本 虽然混合分布的线性函数的后验分布的样本比如上面的预测分布可以直接从折叠采样器的实现中计算出来但是对于非线性函数的推断需要我们首先从混合分布中生成样本。我们可以从随机度量中获得后验样本。需要注意的是为了从 得到后验样本我们需要监控所有参与其计算的随机变量即成员变量xi聚类参数muTilde和s2Tilde以及浓度参数alpha。 点击标题查阅往期内容 WINBUGS对随机波动率模型进行贝叶斯估计与比较 左右滑动查看更多 01 02 03 04 下面的代码从随机测量中生成后验样本。cMCMC对象包括模型和参数的后验样本。函数估计了一个截断水平即truncG。后验样本是一个带列的矩阵其中参数分布向量的维度在本例中为。 outputG - getSamplesDPmeasure(cmcmc) 下面的代码使用随机测量的后验样本来计算的后验样本。请注意这些样本是基于转换后的模型计算的大于70的值对应于上述定义的网格上大于0.035的值。 truncG - outputG$trunc # G的截断水平probY70 - rep(0, nrow(samples)) # P(y.tilde70)的后验样本hist(probY70 ) 使用CRP表示法拟合伽马混合分布 不限于在DPM模型中使用高斯核。就Old Faithful数据而言除了我们在上一节中介绍的对数尺度上的高斯核的混合分布外还有一种选择是数据原始尺度上的伽马混合分布。 模型 在这种情况下模型的形式为 其中对应于两个独立Gamma分布的乘积。下面的代码提供了该模型。 y\[i\] ~ dgamma(shape beta\[i\], scale lambda\[i\])beta\[i\] - betaTilde\[xi\[i\]\]lambda\[i\] - lambdaTilde\[xi\[i\]\] 请注意在这种情况下向量beta和lambda的长度为 。这样做是为了减少与采样算法有关的计算和存储负担。你可以把这种方法看作是对过程的截断只不过它可以被认为是*精确的截断。事实上在CRP表示法下只要采样器的成分数严格低于采样器每次迭代的参数向量的长度使用长度短于样本中观察值的参数向量就会生成一个合适的算法。 运行MCMC算法 下面的代码设置了模型数据和常数初始化了参数定义了模型对象并建立和运行了Gamma混合分布的MCMC算法。请注意在构建MCMC时会产生一个关于聚类参数数量的警告信息。这是因为betaTilde和lambdaTilde的长度小于。另外请注意在执行过程中没有产生错误信息这表明所需的集群数量未超过50个的上限。 data - list(y waiting)
Model(code, data data) cModel - compile samples - runMCMC(cmcmc, niter 7000, nburnin 2000, setSeed TRUE) 在这种情况下我们使用参数的后验样本来构建一个轨迹图并估计的后验分布。 # 参数的后验样本的跟踪图
ts.plot(samples\[ , alpha\], xlab iteration, ylab expression(alpha)) # 参数的后验样本的直方图
hist(samples\[, alpha\]) 从混合分布中生成样本 和以前一样我们从后验分布中获得样本。 outputG - getSamplesDPmeasure(cmcmc) 我们使用这些样本来创建一个数据密度的估计值以及一个95%置信带。 for(iter in seq_len)) {density\[iter, \] - sapply(grid, function(x)sum( weightSamples\[iter, \] * dgamma)))
}hist(waiting, freq FALSE 我们再次看到数据的密度是双峰的看起来与我们之前得到的数据非常相似。 使用stick-breaking 表示法拟合伽马DP混合分布 模型 Dirichlet过程混合物的另一种表示方法是使用随机分布的stick-breaking表示Sethuraman, 1994。 引入辅助变量表明哪个成分产生了每个观测值上一节讨论的Gamma密度的混合物的相应模型的形式为 其中 是两个独立Gamma分布的乘积。 与. 下面的代码提供了该模型说明。 y\[i\] ~ dgamma(shape beta\[i\], scale lambda\[i\])beta\[i\] - betaStar\[z\[i\]\]lambda\[i\] - lambdaStar\[z\[i\]\]z\[i\] ~ dcat(w\[1:Trunc\])# stick-breaking v\[i\] ~ dbeta(1, alpha)w\[1:Trunc\] - stick_breaking(v\[1:(Trunc-1)\]) # stick-breaking 权重betaStar\[i\] ~ dgamma(shape 71, scale 2) 注意截断水平已被设置为Trunc值该值将在函数的常数参数中定义。 运行MCMC算法 下面的代码设置了模型数据和常量初始化了参数定义了模型对象并建立和运行了Gamma混合分布的MCMC算法。当使用stick-breaking表示时会指定一个分块Gibbs抽样器Ishwaran, 2001; Ishwaran and James, 2002。 data - list(y waiting)
consts -length(waiting)
betaStar rgammalambdaStar rgammav rbetaz samplealpha 1 compile(Model) MCMC(rModel, c(w, betaStar, lambdaStar, z, alpha))
comp(mcmc ) MCMC(cmcmc, niter 24000) 使用stick-breaking近似法会自动提供随机分布的近似值即 。下面的代码使用来自样本对象的后验样本计算后验样本并从中计算出数据的密度估计。 densitySamples\[i, \] - sapply(grid, function(x) sum(weightSamples * dgamma(x, shape ,scale )))hist( waiting ylimc(0,0.05), 正如预期的那样这个估计值看起来与我们通过CRP表示的过程获得的估计值相同。 贝叶斯非参数化非参数化随机效应 我们将采用一个参数化的广义线性混合模型并展示如何切换到非参数化的随机效应表示避免了正态分布的随机效应假设。 心肌梗死MIs的参数化meta分析 我们将在对以前非常流行的糖尿病药物 Avandia 的副作用进行meta分析的背景下说明使用非参数混合模型对随机效应分布进行建模。我们分析的数据在引起对这种药物的安全性的严重质疑方面发挥了作用。问题是使用Avandia 是否会增加心肌梗死心脏病发作的风险。每项研究都有治疗和对照组。 模型的制定 我们首先进行基于标准的广义线性混合模型GLMM的meta分析。向量n和x分别包含对照组的患者总数和每项研究中对照组的心肌梗死患者人数。同样向量m和y包含接受药物的病人的类似信息。该模型的形式为 其中随机效应、遵循共同的正态分布、和被合理地赋予非信息性先验。参数量化了对照组和治疗组之间的风险差异而参数则量化了研究的具体变化。 这个模型可以用以下代码指定。 y\[i\] ~ dbin(size m\[i\], prob q\[i\]) # 药物MIsx\[i\] ~ dbin(size n\[i\], prob p\[i\]) # 控制MIsq\[i\] - expit(theta gamma\[i\]) # 药物的对数指数p\[i\] - expit(gamma\[i\]) #对照组对数gamma\[i\] ~ dnorm(mu, var tau2) # 研究效果theta ~ dflat() # 药物的影响# 随机效应超参数mu ~ dnorm(0, 10)tau2 ~ dinvgamma(2, 1) 运行MCMC 让我们来运行一个基本的MCMC。 MCMC(codeParam, data, inits,constants, monitors c(mu, tau2, theta, gamma) par(mfrow c(1, 4)hist(gammaMn)
hist(samples\[1000, gammaCols) 结果表明对照组和治疗组之间存在着整体的风险差异。但是正态性假设呢我们的结论对该假设是否稳健也许随机效应的分布是偏斜的。 用于meta分析的基于DP的随机效应模型 模型 现在我们对使用非参数分布。更具体地说我们假设每个都是由位置尺度的正态混合分布产生的。 这种模型引起了随机效应之间的聚类。与密度估计问题的情况一样DP先验允许数据决定分量的数量从最少的一个分量即简化为参数模型到最多的分量即每个观测值有一个分量。如果数据支持这种行为这允许随机效应的分布是多模态的大大增加了其灵活性。这个模型可以用以下代码指定。 y\[i\] ~ dbin(size m\[i\], prob q\[i\]) # 药物MIsx\[i\] ~ dbin(size n\[i\], prob p\[i\]) # MIsq\[i\] - expit(theta gamma\[i\]) # 药物的对数指数p\[i\] - expit(gamma\[i\]) # 对照组对数值gamma\[i\] ~ dnorm(mu\[i\], var tau2\[i\]) # 来自混合物的随机效应。mu\[i\]- muTilde\[xi\[i\]\] # 来自聚类的随机效应的平均值 xi\[i\]tau2\[i\] - tau2Tilde\[xi\[i\]\] # 来自群组xi\[i\]的随机效应变量# 从基础测量中提取混合成分参数muTilde\[i\] ~ dnorm(mu0, var var0)tau2Tilde\[i\] ~ dinvgamma(a0, b0)# 用于将研究报告聚类为混合成分的CRPxi\[1:nStudies\] ~ dCRP(alpha, size nStudies)# 超参数theta ~ dflat() # 药物的影响 运行MCMC 以下代码对模型进行了编译并对模型运行了一个压缩Gibbs抽样 inits - list(gamma rnorm(nStudies))MCMC(code BNP, data data) hist(samplesBNP\[, theta\], xlab expression(theta), main avandia的影响)main 随机效应分布)main 随机效应分布)# 推断出了多少个混合成分
xiRes - samplesBNP\[, xiCols\]. 主要推论似乎对原始的参数化假设很稳健。这可能是由于没有太多证据表明随机效应分布中缺乏正态性。 参考文献 Blackwell, D. and MacQueen, J. 1973. Ferguson distributions via Polya urn schemes. The Annals of Statistics 1:353-355. Ferguson, T.S. 1974. Prior distribution on the spaces of probability measures. Annals of Statistics 2:615-629. Lo, A.Y. 1984. On a class of Bayesian nonparametric estimates I: Density estimates. The Annals of Statistics 12:351-357. 点击文末“阅读原文” 获取全文完整资料。 本文选自《R语言贝叶斯非参数模型密度估计、非参数化随机效应META分析心肌梗死数据》。 点击标题查阅往期内容 R语言用贝叶斯线性回归、贝叶斯模型平均 (BMA)来预测工人工资 R语言逻辑回归、Naive Bayes贝叶斯、决策树、随机森林算法预测心脏病 PYTHON银行机器学习:回归、随机森林、KNN近邻、决策树、高斯朴素贝叶斯、支持向量机SVM分析营销活动数据|数据分享 PYTHON用户流失数据挖掘建立逻辑回归、XGBOOST、随机森林、决策树、支持向量机、朴素贝叶斯和KMEANS聚类用户画像 MATLAB随机森林优化贝叶斯预测分析汽车燃油经济性 R语言中贝叶斯网络BN、动态贝叶斯网络、线性模型分析错颌畸形数据 使用贝叶斯层次模型进行空间数据分析 MCMC的rstan贝叶斯回归模型和标准线性回归模型比较 python贝叶斯随机过程马尔可夫链Markov-ChainMC和Metropolis-HastingsMH采样算法可视化 Python贝叶斯推断Metropolis-HastingsM-HMCMC采样算法的实现 matlab贝叶斯隐马尔可夫hmm模型实现 贝叶斯线性回归和多元线性回归构建工资预测模型 Metropolis Hastings采样和贝叶斯泊松回归Poisson模型 贝叶斯分位数回归、lasso和自适应lasso贝叶斯分位数回归分析免疫球蛋白、前列腺癌数据 R语言RSTAN MCMCNUTS采样算法用LASSO 构建贝叶斯线性回归模型分析职业声望数据 R语言STAN贝叶斯线性回归模型分析气候变化影响北半球海冰范围和可视化检查模型收敛性 PYTHON用户流失数据挖掘建立逻辑回归、XGBOOST、随机森林、决策树、支持向量机、朴素贝叶斯和KMEANS聚类用户画像 贝叶斯分位数回归、lasso和自适应lasso贝叶斯分位数回归分析免疫球蛋白、前列腺癌数据 R语言JAGS贝叶斯回归模型分析博士生延期毕业完成论文时间 R语言Metropolis Hastings采样和贝叶斯泊松回归Poisson模型 Python决策树、随机森林、朴素贝叶斯、KNNK-最近邻居分类分析银行拉新活动挖掘潜在贷款客户 R语言贝叶斯MCMC用rstan建立线性回归模型分析汽车数据和可视化诊断 R语言贝叶斯MCMCGLM逻辑回归、Rstan线性回归、Metropolis Hastings与Gibbs采样算法实例 R语言贝叶斯Poisson泊松-正态分布模型分析职业足球比赛进球数 随机森林优化贝叶斯预测分析汽车燃油经济性 R语言逻辑回归、Naive Bayes贝叶斯、决策树、随机森林算法预测心脏病 R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数 R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归 Python贝叶斯回归分析住房负担能力数据集 R语言实现贝叶斯分位数回归、lasso和自适应lasso贝叶斯分位数回归分析 Python用PyMC3实现贝叶斯线性回归模型 R语言用WinBUGS 软件对学术能力测验建立层次分层贝叶斯模型 R语言Gibbs抽样的贝叶斯简单线性回归仿真分析 R语言和STAN,JAGS用RSTAN,RJAG建立贝叶斯多元线性回归预测选举数据 R语言基于copula的贝叶斯分层混合模型的诊断准确性研究 R语言贝叶斯线性回归和多元线性回归构建工资预测模型 R语言贝叶斯推断与MCMC实现Metropolis-Hastings 采样算法示例 R语言stan进行基于贝叶斯推断的回归模型 R语言中RStan贝叶斯层次模型分析示例 R语言使用Metropolis-Hastings采样算法自适应贝叶斯估计与可视化 R语言随机搜索变量选择SSVS估计贝叶斯向量自回归BVAR模型 WinBUGS对多元随机波动率模型贝叶斯估计与模型比较 R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样 R语言贝叶斯推断与MCMC实现Metropolis-Hastings 采样算法示例 R语言使用Metropolis-Hastings采样算法自适应贝叶斯估计与可视化 视频R语言中的Stan概率编程MCMC采样的贝叶斯模型 R语言MCMC:Metropolis-Hastings采样用于回归的贝叶斯估计
