原文链接：http://tecdat.cn/?p=24103
此示例说明如何使用逻辑回归模型进行贝叶斯推断 （ 点击文末“阅读原文”获取完整代码数据 ）。
统计推断通常基于最大似然估计 (MLE)。MLE 选择能够使数据似然最大化的参数，是一种较为自然的方法。在 MLE 中，假定参数是未知但固定的数值，并在一定的置信度下进行计算。在贝叶斯统计中，使用概率来量化未知参数的不确定性，因而未知参数被视为随机变量。
贝叶斯推断
贝叶斯推断是结合有关模型或模型参数的先验知识来分析统计模型的过程。这种推断的根基是贝叶斯定理：
例如，假设我们有正态观测值
其中 sigma 是已知的，theta 的先验分布为
在此公式中，mu 和 tau（有时也称为超参数）也是已知的。如果观察 X 的 n 个样本，我们可以获得 theta 的后验分布
下图显示 theta 的先验、似然和后验。
y = norpdf(thta, posMan,psSD);
plot(theta'-', theta,'--', theta,'-.')汽车实验数据
在一些简单的问题中，例如前面的正态均值推断示例，很容易计算出封闭形式的后验分布。但是，在涉及非共轭先验的一般问题中，后验分布很难或不可能通过分析来进行计算。我们将以逻辑回归作为示例。此示例包含一个实验，以帮助建模不同重量的汽车在里程测试中的未通过比例。数据包括被测汽车的重量、汽车数量以及失败次数等观测值。我们采用一组经过变换的重量，以减少回归参数估值的相关性。
% 一组汽车的重量
% 每个重量下测试的汽车数量
[48 42 31 34 31 21 23 23 21 16 17 21]';
% 在每个重量上有不良mpg表现的汽车数量
[1 2 0 3 8 8 14 17 19 15 17 21]';逻辑回归模型
逻辑回归（广义线性模型的一种特例）适合这些数据，因为因变量呈二项分布。逻辑回归模型可以写作：
其中 X 是设计矩阵，b 是包含模型参数的向量。我们可以将此方程写作：
@(b,x) exp(b(1)+b(2).*x)./(1+exp(b(1)+b(2).*x));如果您有一些先验知识或者已经具备某些非信息性先验，则可以指定模型参数的先验概率分布。例如，在此示例中，我们使用正态先验值表示截距 b1 和斜率 b2，即
@(b1) normpdf(b1,0,20); % 截距的先验。
@(b2) normpdf(b2,0,20); % 斜率的先验。根据贝叶斯定理，模型参数的联合后验分布与似然和先验的乘积成正比。
请注意，此模型中后验的归一化常数很难进行分析。但是，即使不知道归一化常数，如果您知道模型参数的大致范围，也可以可视化后验分布。
msh(b2,b1,sipot)
view(-10,30)此后验沿参数空间的对角线伸长，表明（在我们观察数据后）我们认为参数是相关的。这很有意思，因为在我们收集任何数据之前，我们假设它们是独立的。相关性来自我们的先验分布与似然函数的组合。
点击标题查阅往期内容
R语言随机森林RandomForest、逻辑回归Logisitc预测心脏病数据和可视化分析
左右滑动查看更多
01
02
03
04
_切片_采样
蒙特卡罗方法常用于在贝叶斯数据分析中汇总后验分布。其想法是，即使您不能通过分析的方式计算后验分布，也可以从分布中生成随机样本，并使用这些随机值来估计后验分布或推断的统计量，如后验均值、中位数、标准差等。_切片_采样是一种算法，用于从具有任意密度函数的分布中进行抽样，已知项最多只有一个比例常数 - 而这正是从归一化常数未知的复杂后验分布中抽样所需要的。此算法不生成独立样本，而是生成马尔可夫序列，其平稳分布就是目标分布。因此，切片抽样器是一种马尔可夫链蒙特卡罗 (MCMC) 算法。但是，它与其他众所周知的 MCMC 算法不同，因为只需要指定缩放的后验，不需要建议分布或边缘分布。
此示例说明如何使用切片抽样器作为里程测试逻辑回归模型的贝叶斯分析的一部分，包括从模型参数的后验分布生成随机样本、分析抽样器的输出，以及对模型参数进行推断。第一步是生成随机样本。
sliesmle(inial,nsapes,'pdf');采样器输出分析
从切片采样获取随机样本后，很重要的一点是研究诸如收敛和混合之类的问题，以确定将样本视为是来自目标后验分布的一组随机实现是否合理。观察边缘轨迹图是检查输出的最简单方法。
plot(trace(:,1))从这些图中可以明显看出，在处理过程趋于平稳之前，参数起始值的影响会维持一段时间（大约 50 个样本）才会消失。
检查收敛以使用移动窗口计算统计量（例如样本的均值、中位数或标准差）也很有帮助。这样可以产生比原始样本轨迹更平滑的图，并且更容易识别和理解任何非平稳性。
mvag = fier( (1/50)*os(50,1), 1, tace);
plot(moav(:,1))由于这些是基于包含 50 次迭代的窗口计算的移动平均值，因此前 50 个值无法与图中的其他值进行比较。然而，每个图的其他值似乎证实参数后验均值在 100 次左右迭代后收敛至平稳分布。同样显而易见的是，这两个参数彼此相关，与之前的后验密度图一致。
由于磨合期代表目标分布中不能合理视为随机实现的样本，因此不建议使用切片采样器一开始输出的前 50 个左右的值。您可以简单地删除这些输出行，但也可以指定一个“预热”期。在已知合适的预热长度（可能来自先前的运行）时，这种方式很简便。
slcsapl(inial,nsmes,'pf',pot, ..'brin',50);
plot(trace(:,1))这些跟踪图没有显示出任何不平稳，表明预热期已完成。
但是，还需要了解跟踪图的另一方面。虽然截距的轨迹看起来像高频噪声，但斜率的轨迹好像具有低频分量，表明相邻迭代的值之间存在自相关。虽然也可以从这个自相关样本计算均值，但我们通常会通过删除样本中的冗余数据这一简便的操作来降低存储要求。如果它同时消除了自相关，我们还可以将这些数据视为独立值样本。例如，您可以通过只保留第 10 个、第 20 个、第 30 个等值来稀释样本。
sceampe(...
'brin'50,'tin',10);要检查这种稀释的效果，可以根据轨迹估计样本自相关函数，并使用它们来检查样本是否快速混合。
fftetendtrce,'cnsant');
F .* coj(F);
for i = 1:2
lineles = stem(:20, F(:i) 'filled' , 'o');第一个滞后的自相关值对于截距参数很明显，对于斜率参数更是如此。我们可以使用更大的稀释参数重复抽样，以进一步降低相关性。但为了完成本示例的目的，我们将继续使用当前样本。
推断模型参数
与预期相符，样本直方图模拟了后验密度图。
hist(rce,[25,25]);
view(-10,30)您可以使用直方图或核平滑密度估计值来总结后验样本的边缘分布属性。
kdeiy(rae(:2))您还可以计算描述性统计量，例如随机样本的后验均值或百分位数。为了确定样本大小是否足以实现所需的精度，将所需的轨迹统计量作为样本数的函数来进行查看会很有帮助。
csu= csm(rae);
plot(csm(:,1)'./(1:sals))在这种情况下，样本大小 1000 似乎足以为后验均值估计值提供良好的精度。
mean(te)总结
您能够轻松地指定似然和先验。您也可以将它们结合起来用于推断后验分布。您可以通过马尔可夫链蒙特卡罗仿真在 MATLAB 中执行贝叶斯分析。
点击文末 “阅读原文”
获取全文完整资料。
本文选自《matlab用马尔可夫链蒙特卡罗 (MCMC) 的Logistic逻辑回归模型分析汽车实验数据》。
点击标题查阅往期内容
R语言coda贝叶斯MCMC Metropolis-Hastings采样链分析和收敛诊断可视化
R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样
R语言贝叶斯METROPOLIS-HASTINGS GIBBS 吉布斯采样器估计变点指数分布分析泊松过程车站等待时间
R语言马尔可夫MCMC中的METROPOLIS HASTINGS，MH算法抽样（采样）法可视化实例
python贝叶斯随机过程：马尔可夫链Markov-Chain，MC和Metropolis-Hastings，MH采样算法可视化
Python贝叶斯推断Metropolis-Hastings（M-H）MCMC采样算法的实现
Metropolis Hastings采样和贝叶斯泊松回归Poisson模型
Matlab用BUGS马尔可夫区制转换Markov switching随机波动率模型、序列蒙特卡罗SMC、M H采样分析时间序列R语言RSTAN MCMC：NUTS采样算法用LASSO 构建贝叶斯线性回归模型分析职业声望数据
R语言BUGS序列蒙特卡罗SMC、马尔可夫转换随机波动率SV模型、粒子滤波、Metropolis Hasting采样时间序列分析
R语言Metropolis Hastings采样和贝叶斯泊松回归Poisson模型
R语言贝叶斯MCMC：用rstan建立线性回归模型分析汽车数据和可视化诊断
R语言贝叶斯MCMC：GLM逻辑回归、Rstan线性回归、Metropolis Hastings与Gibbs采样算法实例
R语言贝叶斯Poisson泊松-正态分布模型分析职业足球比赛进球数
R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数
R语言逻辑回归、Naive Bayes贝叶斯、决策树、随机森林算法预测心脏病
R语言中贝叶斯网络（BN）、动态贝叶斯网络、线性模型分析错颌畸形数据
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采样用于回归的贝叶斯估计R语言逻辑回归、Naive Bayes贝叶斯、决策树、随机森林算法预测心脏病
R语言逻辑回归(Logistic Regression)、回归决策树、随机森林信用卡违约分析信贷数据集PYTHON用户流失数据挖掘：建立逻辑回归、XGBOOST、随机森林、决策树、支持向量机、朴素贝叶斯和KMEANS聚类用户画像
Python对商店数据进行lstm和xgboost销售量时间序列建模预测分析
PYTHON集成机器学习：用ADABOOST、决策树、逻辑回归集成模型分类和回归和网格搜索超参数优化
R语言集成模型：提升树boosting、随机森林、约束最小二乘法加权平均模型融合分析时间序列数据
Python对商店数据进行lstm和xgboost销售量时间序列建模预测分析
R语言用主成分PCA、 逻辑回归、决策树、随机森林分析心脏病数据并高维可视化
R语言基于树的方法：决策树，随机森林，Bagging，增强树
R语言用逻辑回归、决策树和随机森林对信贷数据集进行分类预测
spss modeler用决策树神经网络预测ST的股票
R语言中使用线性模型、回归决策树自动组合特征因子水平
R语言中自编基尼系数的CART回归决策树的实现
R语言用rle，svm和rpart决策树进行时间序列预测
python在Scikit-learn中用决策树和随机森林预测NBA获胜者
python中使用scikit-learn和pandas决策树进行iris鸢尾花数据分类建模和交叉验证
R语言里的非线性模型：多项式回归、局部样条、平滑样条、 广义相加模型GAM分析
R语言用标准最小二乘OLS，广义相加模型GAM ，样条函数进行逻辑回归LOGISTIC分类
R语言ISLR工资数据进行多项式回归和样条回归分析
R语言中的多项式回归、局部回归、核平滑和平滑样条回归模型
R语言用泊松Poisson回归、GAM样条曲线模型预测骑自行车者的数量
R语言分位数回归、GAM样条曲线、指数平滑和SARIMA对电力负荷时间序列预测R语言样条曲线、决策树、Adaboost、梯度提升(GBM)算法进行回归、分类和动态可视化
如何用R语言在机器学习中建立集成模型？
R语言ARMA-EGARCH模型、集成预测算法对SPX实际波动率进行预测在python 深度学习Keras中计算神经网络集成模型R语言ARIMA集成模型预测时间序列分析R语言基于Bagging分类的逻辑回归(Logistic Regression)、决策树、森林分析心脏病患者
R语言基于树的方法：决策树，随机森林，Bagging，增强树
R语言基于Bootstrap的线性回归预测置信区间估计方法
R语言使用bootstrap和增量法计算广义线性模型（GLM）预测置信区间
R语言样条曲线、决策树、Adaboost、梯度提升(GBM)算法进行回归、分类和动态可视化
Python对商店数据进行lstm和xgboost销售量时间序列建模预测分析
R语言随机森林RandomForest、逻辑回归Logisitc预测心脏病数据和可视化分析
R语言用主成分PCA、 逻辑回归、决策树、随机森林分析心脏病数据并高维可视化
Matlab建立SVM，KNN和朴素贝叶斯模型分类绘制ROC曲线
matlab使用分位数随机森林（QRF）回归树检测异常值