目标

2022/4/17-2022/5/10

实现自适应的MCMC方法（Adaptive Metropolis Algorithm）

本地目录：E:\Research\OptA\MCMC

如有问题，欢迎交流探讨！ 邮箱：lujiabo@hhu.edu.cn 卢家波
来信请说明博客标题及链接，谢谢。

MCMC简介

MCMC方法是基于贝叶斯理论框架，通过建立平衡分布为π(x)\pi(x)π(x)的马尔可夫链，并对其平衡分布进行采样，通过不断更新样本信息而使马尔可夫链能充分搜索模型参数空间，最终收敛于高概率密度区，因此，MCMC方法是对理想的贝叶斯推断过程的一种近似。MCMC方法的关键是如何构造有效的推荐分布，确保按照推荐分布抽取的样本收敛于高概率密度区。

MCMC方法论为建立实际的统计模型提供了一种有效工具，能将一些复杂的高维问题转化为一系列简单的低维问题，因此适用于复杂统计模型的贝叶斯计算。目前，在贝叶斯分析中应用最为广泛的MCMC方法主要基于Gibbs采样方法和Metropolis-Hastings 方法[2,8-9]。在此基础上，后人对这些方法开展了进一步的改进和推广研究，Haario等人(2001)[1]提出了自适应的MCMC方法（Adaptive Metropolis Algorithm），其主要特点是将推荐分布定义为参数空间多维正态分布，在进化过程中自适应地调整协方差矩阵，从而大大提高算法的收敛速度。

模型描述

线性模型

利用AM-MCMC估计线性回归参数，模型如下：

y=kx+by=kx+by=kx+b

参数

进行AM-MCMC不确定性估计的2个参数为 kkk 和 bbb

参数名	含义	真值	下限	上限
kkk	斜率	2	-3	3
bbb	截距	-1	-3	3

实测值

实测值取-5到5之间的整数对应的函数值，每个点增加正态随机数N(0, 0.16)扰动，即

y=−11,−9,−7,−5,−3,−1,1,3,5,7,9,x∈[−5,5]→y=-11,-9,-7,-5, -3, -1, 1, 3, 5, 7, 9, x\in[-5,5]\rightarrowy=−11,−9,−7,−5,−3,−1,1,3,5,7,9,x∈[−5,5]→

y=−10.79,−9.4,−6.77,−5.3,−3.52,−0.97,1.3,3.13,5.16,7.56,8.8,x∈[−5,5]y=-10.79, -9.4, -6.77, -5.3, -3.52, -0.97, 1.3, 3.13, 5.16, 7.56, 8.8,x\in[-5,5]y=−10.79,−9.4,−6.77,−5.3,−3.52,−0.97,1.3,3.13,5.16,7.56,8.8,x∈[−5,5]

AM-MCMC参数设置

参数初始协方差矩阵C0C_0C0为对角矩阵，方差取参数取值范围的1/201/201/20；初始迭代次数t0=100t_0=100t0=100。算法平行运行5次，每次采样5000个样本。

模拟次数	参数个数	并行抽样次数	置信水平%	热身期
5000	2	5	90	500

似然函数

采用BOX&Tiao给出的似然函数，参考[5]p.64(4-14)

p(θt∣y)∝[∑i=1Ne(θt)i2]−12Np(\theta^{t}|y) \propto [\sum_{i=1}^{N}e(\theta^{t})_i^2]^{-\frac{1}{2}N}p(θt∣y)∝[i=1∑Ne(θt)i2]−21N

分析结果

收敛判断

针对单序列是否稳定，可采用平均值法和方差法，考察迭代过程中的参数平均值和方差是否稳定[3]。针对多序列是否收敛，采用Gelman[7]在1992年提出的收敛诊断指标R\sqrt{R}R——比例缩小得分，计算每一参数的比例缩小得分R\sqrt{R}R，若接近于1表示参数收敛到了后验分布上。Gelman和Rubin建议将R<1.2\sqrt{R}<1.2R<1.2作为多序列抽样算法收敛判断条件。

单序列

单一采样序列的平均值和方差的变化过程图。从图可看出，当t>500t>500t>500后，参数kkk和bbb的平均值和方差基本达到稳定，因此单一序列是收敛的。

斜率k

截距b

多序列

R\sqrt{R}R变化过程，在迭代初期，R\sqrt{R}R变化剧烈；当t>200t>200t>200后，R\sqrt{R}R趋于稳定，并接近于1.0，说明参数kkk和bbb的MCMC采样序列均能稳定收敛到其参数的后验分布，且算法全局收敛。综合考虑上述单序列和多序列参数的收敛情况，将每一序列的前500次舍去，这样５次平行试验共采集了22500个样本用于参数后验分布的统计分析。

斜率k

截距b

参数后验分布

根据上述MCMC抽样得到参数kkk和bbb的后验分布。由图可知，bbb比kkk具有更大的不确定性。

斜率k

截距b

置信区间

将似然函数归一化权重根据抽样结果的升序排列，计算5%、50%、95%置信度所在模拟结果，绘制如下图。大多数实测数据都包含于90%的置信区间内。

结论

与Metropolis-Hastings算法相比，Adaptive Metropolis(AM)算法不需要事先确定推荐分布；与Gibbs采样相比，AM不需要导出条件分布，对于复杂的水文模型，通常很难导出单个参数的条件分布；与GLUE方法相比，AM基于贝叶斯理论框架，是对理想贝叶斯推断过程的一种近似，而GLUE主要通过随机抽样估计参数不确定性。

参考文献

[1]Haario, H., Saksman, E., & Tamminen, J. (2001). An Adaptive Metropolis Algorithm. Bernoulli, 7(2), 223–242. https://doi.org/10.2307/3318737
[2]Kuczera G, Parent E. (1998). Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the Metropolis algorithm. Journal of Hydrology, 211(1), 69-85. https://doi.org/10.1016/S0022-1694(98)00198-X
[3]梁忠民,戴荣,綦晶. 基于MCMC的水文模型参数不确定性及其对预报的影响分析[C]. 环境变化与水安全——第五届中国水论坛论文集.,2007:126-130.
[4]李向阳. 水文模型参数优选及不确定性分析方法研究[D].大连理工大学,2006.
[5]曹飞凤. 基于MCMC方法的概念性流域水文模型参数优选及不确定性研究[D].浙江大学,2010.
[6]Thiemann M, Trosset M, Gupta H, Sorooshian S. (2001). Bayesian recursive parameter estimation for hydrologic models. Water Resources Research, 37(10), 2521-2535. https://doi.org/10.1029/2000WR900405
[7]Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical science, 7(4), 457-472. https://doi.org/10.1214/ss/1177011136
[8]W. K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika, Volume 57, Issue 1, April 1970, Pages 97–109, https://doi.org/10.1093/biomet/57.1.97
[9]MetropoUs, N., Rosenbluth, A., Rosenbluth, M., Teller, A., & Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chem. Physics, 21, 1087-1092.https://doi.org/10.1063/1.1699114
[10]Armadillo. C++ library for linear algebra & scientific computing. http://arma.sourceforge.net/
[11]Conrad Sanderson and Ryan Curtin. Armadillo: a template-based C++ library for linear algebra. Journal of Open Source Software, Vol. 1, pp. 26, 2016. http://dx.doi.org/10.21105/joss.00026
[12]Conrad Sanderson and Ryan Curtin. An Adaptive Solver for Systems of Linear Equations. International Conference on Signal Processing and Communication Systems, pp. 1-6, 2020.
[13]Matplotlib for C++. https://github.com/lava/matplotlib-cpp

源码

AM为单次抽样程序，PAM为平行抽样程序，继承于AM类。由于高度耦合，因此AM类成员都设置为公开public。main.cpp为主函数，负责创建PAM类，调用参数估计Estimate()函数。参数的上下限、初始值、实测值、算法参数在AM类的Initialize()函数中设置。平行抽样、置信度在PAM类的构造函数中设置。

程序中用到了Matplotlib-cpp[13]，用于数据可视化，使用方法参见【C++】11 Visual Studio 2019 C++安装matplotlib-cpp绘图。

还用到了线性代数库Armadillo[10-12]，用于多元正态分布抽样，使用方法参见【C++】13 多元正态分布抽样。

AM.h

#pragma once/******************************************************************************
文件名: AM.h
作者: 卢家波
单位：河海大学水文水资源学院
邮箱：lujiabo@hhu.edu.cn
QQ：1847096852
版本：2022.5 创建 V1.0
版权: MIT
引用格式：卢家波，AM-MCMC算法C++实现. 南京：河海大学，2022.LU Jiabo, Adaptive Metropolis algorithm of Markov Chain Monte Carlo in C++. Nanjing:Hohai University, 2022.
参考文献：[1]Haario, H., Saksman, E., & Tamminen, J. (2001). An Adaptive Metropolis Algorithm. Bernoulli, 7(2), 223–242. https://doi.org/10.2307/3318737[2]Kuczera G, Parent E. (1998). Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the Metropolis algorithm. Journal of Hydrology, 211(1), 69-85. https://doi.org/10.1016/S0022-1694(98)00198-X[3]梁忠民,戴荣,綦晶. 基于MCMC的水文模型参数不确定性及其对预报的影响分析[C]. 环境变化与水安全——第五届中国水论坛论文集.,2007:126-130.[4]李向阳. 水文模型参数优选及不确定性分析方法研究[D].大连理工大学,2006.[5]曹飞凤. 基于MCMC方法的概念性流域水文模型参数优选及不确定性研究[D].浙江大学,2010.[6]Thiemann M, Trosset M, Gupta H, Sorooshian S. (2001). Bayesian recursive parameter estimation for hydrologic models. Water Resources Research, 37(10), 2521-2535.  https://doi.org/10.1029/2000WR900405[7]Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical science, 7(4), 457-472. https://doi.org/10.1214/ss/1177011136[8]W. K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika, Volume 57, Issue 1, April 1970, Pages 97–109, https://doi.org/10.1093/biomet/57.1.97[9]MetropoUs, N., Rosenbluth, A., Rosenbluth, M., Teller, A., & Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chem. Physics, 21, 1087-1092.https://doi.org/10.1063/1.1699114[10]Armadillo. C++ library for linear algebra & scientific computing. http://arma.sourceforge.net/[11]Conrad Sanderson and Ryan Curtin. Armadillo: a template-based C++ library for linear algebra. Journal of Open Source Software, Vol. 1, pp. 26, 2016. http://dx.doi.org/10.21105/joss.00026[12]Conrad Sanderson and Ryan Curtin. An Adaptive Solver for Systems of Linear Equations. International Conference on Signal Processing and Communication Systems, pp. 1-6, 2020.[13]Matplotlib for C++. https://github.com/lava/matplotlib-cpp******************************************************************************/#include <vector>
#include <armadillo>   //线性代数运算库class AM
{public://AM-MCMC不确定性估计函数virtual void Estimate();//===============1.输入估计参数及实测值===============void Setup();//====================2.运行模型====================void Run();//private://===================3.不确定性估计===================void Analysis();//===================4.输出估计结果===================void Save();//1.初始化，t=0void Initialize();//2.计算协方差矩阵Ctvoid CalculateCovariance();//3.从推荐分布抽样获得候选点void Sample();//4.计算接受概率alphavoid CalAcceptanceProbability();//5.产生[0, 1]随机数uvoid GenerateRandomNumber();//6.判断是否接受候选点void IfAcceptCandidatePoint();//7.重复2-6直到产生足够的样本为止void Repeat();//计算经验协方差矩阵void CalEmpiricalCovariance();//协方差递推公式void RecurseCovariance();//候选点参数检查bool CheckBound();//目标概率密度函数,取对数double logPDF();//调用待分析模型，得到模拟值void CallModel();//计算似然函数值void CalculateLikeliFun();//计算纳什系数void CalculateNSE();//计算均方误差void CalculateMSE();//计算BOX&Tiao给出的似然函数，参考[5]p.64(4-14)void CalculateBOX();//收敛判断指标void CheckConvergence();//参数后验分布频率直方图virtual void PlotPostDistribution();//绘制置信区间virtual void PlotConfidenceInterval();int dimension;   //参数的空间维度int t0;          //初始迭代次数int iterationNumber;  //迭代总次数int currentIterNum;   //当前迭代次数int burnInNumber;     //热身次数double sd;       //比例因子，仅取决于参数的空间维度double epsilon;  //为一个较小的数，以确保Ci不为奇异矩阵double alpha;    //候选点接受概率double logDensity;   //目标分布的对数密度，$\pi(x)$double candidatePointLogDensity;  //候选点的对数密度double randomNumber;  //随机数//实测值std::vector<double> observedData;//每次待分析模型运行得到的模拟值、模拟值数组std::vector<double> modelledData;std::vector< std::vector<double> > modelledDataArray;//似然函数值std::vector<double> NSE;   //纳什效率系数std::vector<double> MSE;   //均方误差std::vector<double> BOX;   //BOX似然函数std::vector<std::string> paramName;  //参数名arma::vec lowerBound;   //参数下限arma::vec upperBound;   //参数上限arma::vec candidatePoint;  //候选点arma::vec currentMean;   //前t次迭代样本点的均值arma::vec previousMean;  //前t-1次迭代样本点的均值arma::vec sample;     //当前样本点arma::mat sampleGroup;  //样本点集合arma::mat sampleMean;   //样本点各次迭代均值arma::mat sampleVariance;  //样本点各次迭代方差arma::mat sampleGroupBurnIn;  //除去热身期的样本arma::mat C0;    //初始协方差矩阵arma::mat Ct;    //第t次迭代协方差矩阵arma::mat Id;    //单位矩阵arma::mat empiricalCovariance;  //cov(X0, ..., Xt0)，经验协方差矩阵};

PAM.h

#pragma once//Parallel Adaptive Metropolis algorithm#include "AM.h"
#include <armadillo>   //线性代数运算库class PAM : public AM
{public:void Estimate() override;PAM();private:void ParallelRun();void CalculateSRF();void PlotSRF();void PlotPostDistribution() override;void PlotConfidenceInterval() override;void Plot();//计算似然函数值对应的归一化权重void CalculateWeight();//经验累积分布函数（CDF）的预测void ecdfPred();//按模拟值的升序对一组权重数组和对应的模拟值数组进行排序void sort(std::vector<double>& w, std::vector<double>& x);//计算权重累加之和std::vector<double> cumsum(const std::vector<double>& w);//计算各百分位数下的样本索引值std::vector<double> CalPercentiles(const std::vector<double>& modValue);//通过置信水平计算百分位数void SetPercentile();int parallelSampleNumber;  //并行抽样次数int modelledDataSize;  //模拟值数据大小double confidenceLevel; //置信水平std::vector<double> parallelBOXBurnIn;  //并行抽样除去热身期的BOX似然函数std::vector<double> weight;std::vector< std::vector<double> > modelledDataArrayBurnIn;  //去除热身期模拟值数组//不确定性边界预测std::vector<double> percentile;  //置信区间的百分位数std::vector<double> ecdf;        //经验累积分布函数std::vector<double> sampLowerPtile;  //低分位数std::vector<double> sampUpperPtile;  //高分位数std::vector<double> sampMedian;      //50%分位数arma::mat scaleReductionFactor;  //比例缩小得分SRF$\sqrt{R}$arma::cube parallelSampleMean;   //并行抽样样本点各次迭代均值arma::cube parallelSampleVariance;  //并行抽样各次迭代方差arma::cube parallelSampleGroupBurnIn;  //并行抽样除去热身期的样本};

AM.cpp

/******************************************************************************
文件名: AM.cpp
作者: 卢家波
单位：河海大学水文水资源学院
邮箱：lujiabo@hhu.edu.cn
QQ：1847096852
版本：2022.5 创建 V1.0
版权: MIT
引用格式：卢家波，AM-MCMC算法C++实现. 南京：河海大学，2022.LU Jiabo, Adaptive Metropolis algorithm of Markov Chain Monte Carlo in C++. Nanjing:Hohai University, 2022.
参考文献：[1]Haario, H., Saksman, E., & Tamminen, J. (2001). An Adaptive Metropolis Algorithm. Bernoulli, 7(2), 223–242. https://doi.org/10.2307/3318737[2]Kuczera G, Parent E. (1998). Monte Carlo assessment of parameter uncertainty in conceptual catchment models: the Metropolis algorithm. Journal of Hydrology, 211(1), 69-85. https://doi.org/10.1016/S0022-1694(98)00198-X[3]梁忠民,戴荣,綦晶. 基于MCMC的水文模型参数不确定性及其对预报的影响分析[C]. 环境变化与水安全——第五届中国水论坛论文集.,2007:126-130.[4]李向阳. 水文模型参数优选及不确定性分析方法研究[D].大连理工大学,2006.[5]曹飞凤. 基于MCMC方法的概念性流域水文模型参数优选及不确定性研究[D].浙江大学,2010.[6]Thiemann M, Trosset M, Gupta H, Sorooshian S. (2001). Bayesian recursive parameter estimation for hydrologic models. Water Resources Research, 37(10), 2521-2535.  https://doi.org/10.1029/2000WR900405[7]Gelman, A., & Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical science, 7(4), 457-472. https://doi.org/10.1214/ss/1177011136[8]W. K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika, Volume 57, Issue 1, April 1970, Pages 97–109, https://doi.org/10.1093/biomet/57.1.97[9]MetropoUs, N., Rosenbluth, A., Rosenbluth, M., Teller, A., & Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chem. Physics, 21, 1087-1092.https://doi.org/10.1063/1.1699114[10]Armadillo. C++ library for linear algebra & scientific computing. http://arma.sourceforge.net/[11]Conrad Sanderson and Ryan Curtin. Armadillo: a template-based C++ library for linear algebra. Journal of Open Source Software, Vol. 1, pp. 26, 2016. http://dx.doi.org/10.21105/joss.00026[12]Conrad Sanderson and Ryan Curtin. An Adaptive Solver for Systems of Linear Equations. International Conference on Signal Processing and Communication Systems, pp. 1-6, 2020.[13]Matplotlib for C++. https://github.com/lava/matplotlib-cpp******************************************************************************/#include "AM.h"#include <random>
#include <numeric>
#include <armadillo>   //线性代数运算库
#include "matplotlibcpp.h"  //matplotlib的C++接口，用法见【C++】11 Visual Studio 2019 C++安装matplotlib-cpp绘图（https://blog.csdn.net/weixin_43012724/article/details/124051588）
namespace plt = matplotlibcpp;double AM::logPDF()
{CallModel();CalculateLikeliFun();//BOX效果最好，MSE次之，NSE最差//if (NSE.back() <= 0)//{//    return log(1.0e-6);//}//else//{//   return log(NSE.back());//}//return log(1.0 / MSE.back());return BOX.back();}void AM::CallModel()
{//测试模型为线性模型 y = k * x + b, x in [-5, 5]//k 真值为 2，b 真值为 -1for (int x = -5; x <= 5; ++x){double temp = candidatePoint(0) * x + candidatePoint(1);modelledData.emplace_back(temp);}modelledDataArray.emplace_back(modelledData);
}void AM::CalculateLikeliFun()
{CalculateNSE();CalculateMSE();CalculateBOX();modelledData.clear();
}void AM::CalculateNSE()
{double measuredValuesSum = std::accumulate(observedData.begin(), observedData.end(), 0.0);double measuredValuesAvg = measuredValuesSum / observedData.size();double numerator = 0.0;double denominator = 0.0;for (double val : observedData){denominator += pow(val - measuredValuesAvg, 2);}for (int index = 0; index < observedData.size(); ++index){numerator += pow(modelledData.at(index) - observedData.at(index), 2);}double nse = 1 - numerator / denominator;NSE.emplace_back(nse);
}void AM::CalculateMSE()
{double squareErrorSum = 0;int dataSize = observedData.size();for (int i = 0; i < dataSize; ++i){squareErrorSum += pow(modelledData.at(i) - observedData.at(i), 2);}double mse = squareErrorSum / static_cast<double>(dataSize);MSE.emplace_back(mse);
}void AM::CalculateBOX()
{double squareErrorSum = 0;int dataSize = observedData.size();for (int i = 0; i < dataSize; ++i){squareErrorSum += pow(modelledData.at(i) - observedData.at(i), 2);}double box = log(pow(squareErrorSum, -static_cast<double>(dataSize) / 2.0));BOX.emplace_back(box);}void AM::CheckConvergence()
{for (int i = 0; i < dimension; ++i){arma::rowvec paramMeanVec = sampleMean.row(i);  //参数均值std::vector<double> paramMean //将arma::vec转化为std::vector<double>类型= arma::conv_to< std::vector<double> >::from(paramMeanVec);  plt::plot(paramMean);plt::title("Change process diagram of the mean value of parameter " + paramName.at(i));plt::xlabel("Number of simulations");plt::ylabel("Mean of " + paramName.at(i));plt::save(paramName.at(i) + "_Mean" + ".png");plt::show();arma::rowvec paramVarianceVec = sampleVariance.row(i);   //参数方差std::vector<double> paramVariance = arma::conv_to< std::vector<double> >::from(paramVarianceVec);plt::plot(paramVariance);plt::title("Change process diagram of the variance value of parameter " + paramName.at(i));plt::xlabel("Number of simulations");plt::ylabel("Variance of " + paramName.at(i));plt::save(paramName.at(i) + "_Variance" + ".png");plt::show();}}void AM::PlotPostDistribution()
{sampleGroupBurnIn = sampleGroup.cols(burnInNumber, iterationNumber);for (int i = 0; i < dimension; ++i){arma::rowvec paramChainVec = sampleGroupBurnIn.row(i);std::vector<double> paramChain = arma::conv_to< std::vector<double> >::from(paramChainVec);plt::hist(paramChain);plt::show();}}void AM::PlotConfidenceInterval()
{}void AM::Estimate()
{Setup();Run();Analysis();Save();
}void AM::Setup()
{Initialize();
}void AM::Run()
{Repeat();
}void AM::Analysis()
{CheckConvergence();PlotPostDistribution();PlotConfidenceInterval();
}void AM::Save()
{sampleGroup.save("sampleGroup.csv", arma::csv_ascii);sampleGroupBurnIn.save("sampleGroupBurbIn.csv", arma::csv_ascii);sampleMean.save("sampleMean.csv", arma::csv_ascii); sampleVariance.save("sampleVariance.csv", arma::csv_ascii);std::ofstream fout("likelihood.csv");fout << "NSE" << "," << "MSE" << "," << "logBOX" << ",\n";int likelihoodSize = BOX.size();for (int i = 0; i < likelihoodSize; ++i){fout << NSE.at(i) << "," << MSE.at(i) << "," << BOX.at(i) << ",\n";}fout.close();
}void AM::Initialize()
{dimension = 2;t0 = 100;iterationNumber = 5000;currentIterNum = 0;burnInNumber = 500;sd = 2.4 * 2.4 / dimension;epsilon = 1.0e-5;alpha = 0;candidatePointLogDensity = 0;randomNumber = 0;observedData = { -10.79, -9.4, -6.77, -5.3, -3.52, -0.97, 1.3, 3.13, 5.16, 7.56, 8.8 };  //每个点增加正态随机数N(0, 0.16)扰动paramName = { "k", "b" };lowerBound = {-3, -3};upperBound = {3, 3 };candidatePoint.zeros(dimension);currentMean = { 1, 1 };previousMean.zeros(dimension);sample = {1, 1};logDensity = logPDF();sampleGroup.zeros(dimension, iterationNumber + 1);sampleGroup.col(0) = sample;sampleMean.zeros(dimension, iterationNumber + 1);sampleMean.col(0) = currentMean;sampleVariance.zeros(dimension, iterationNumber + 1);C0.zeros(dimension, dimension);C0.diag() = (upperBound - lowerBound) / 20.0;Ct = C0;Id.eye(dimension, dimension);empiricalCovariance.zeros(dimension, dimension);
}void AM::CalculateCovariance()
{if (currentIterNum == t0 + 1){CalEmpiricalCovariance();Ct = sd * empiricalCovariance + sd * epsilon * Id;}if (currentIterNum > t0 + 1){RecurseCovariance();}
}void AM::Sample()
{candidatePoint = arma::mvnrnd(sample, Ct);
}void AM::CalAcceptanceProbability()
{if (CheckBound()){candidatePointLogDensity = logPDF();double temp = candidatePointLogDensity - logDensity;alpha = std::min(1.0, exp(temp));}else{alpha = 0.0;modelledDataArray.emplace_back(modelledDataArray.back());NSE.emplace_back(NSE.back());MSE.emplace_back(MSE.back());BOX.emplace_back(BOX.back());}}void AM::GenerateRandomNumber()
{std::random_device rd;std::mt19937 gen(rd());std::uniform_real_distribution<double> u(0.0, 1.0);randomNumber = u(gen);
}void AM::IfAcceptCandidatePoint()
{if (alpha >= randomNumber){sample = candidatePoint;logDensity = candidatePointLogDensity;}sampleGroup.col(currentIterNum) = sample;//更新方差，注意不要对整个sampleGroup求方差，因为包含初始化的0sampleVariance.col(currentIterNum) = arma::var(sampleGroup.cols(0, currentIterNum), 0, 1);  //对矩阵每一行求方差, N-1//更新均值previousMean = currentMean;currentMean = (currentIterNum * previousMean + sample) / static_cast<double>(currentIterNum + 1);sampleMean.col(currentIterNum) = currentMean;
}void AM::Repeat()
{while (currentIterNum < iterationNumber){currentIterNum++;CalculateCovariance();Sample();CalAcceptanceProbability();GenerateRandomNumber();IfAcceptCandidatePoint();}sampleGroupBurnIn = sampleGroup.cols(burnInNumber, iterationNumber);
}void AM::CalEmpiricalCovariance()
{arma::mat tempSum(dimension, dimension, arma::fill::zeros);for (int col = 0; col < sampleGroup.n_cols; col ++){arma::vec x = sampleGroup.col(col);tempSum += x * x.t();}empiricalCovariance = (tempSum - (t0 + 1) * previousMean * previousMean.t())/ static_cast<double>(t0);
}void AM::RecurseCovariance()
{int t = currentIterNum;Ct = (static_cast<double>(t - 1) / static_cast<double>(t)) * Ct+ sd / static_cast<double>(t)* (t * previousMean * previousMean.t()- (t + 1) * currentMean * currentMean.t()+ sample * sample.t()+ epsilon * Id);
}bool AM::CheckBound()
{for (size_t i = 0; i < dimension; i++){if ((candidatePoint(i) < lowerBound(i)) || (candidatePoint(i) > upperBound(i))){return false;}}return true;
}

PAM.cpp

#include "PAM.h"#include <algorithm>
#include "matplotlibcpp.h"  //matplotlib的C++接口，用法见【C++】11 Visual Studio 2019 C++安装matplotlib-cpp绘图（https://blog.csdn.net/weixin_43012724/article/details/124051588）
namespace plt = matplotlibcpp;void PAM::Estimate()
{ParallelRun();CheckConvergence();CalculateSRF();PlotSRF();PlotPostDistribution();PlotConfidenceInterval();Save();
}PAM::PAM()
{Setup();confidenceLevel = 0.9;parallelSampleNumber = 5;scaleReductionFactor = arma::mat(dimension, iterationNumber, arma::fill::zeros);  //比例缩小得分SRF$\sqrt{R}$parallelSampleMean = arma::cube(dimension, iterationNumber + 1, parallelSampleNumber);   //并行抽样样本点各次迭代均值parallelSampleVariance = arma::cube(dimension, iterationNumber + 1, parallelSampleNumber);  //并行抽样各次迭代方差parallelSampleGroupBurnIn = arma::cube(dimension, iterationNumber - burnInNumber + 1, parallelSampleNumber);  //并行抽样除去热身期的样本
}void PAM::ParallelRun()
{for (int i = 0; i < parallelSampleNumber; ++i){Setup();Run();parallelSampleMean.slice(i) = sampleMean;parallelSampleVariance.slice(i) = sampleVariance;parallelSampleGroupBurnIn.slice(i) = sampleGroupBurnIn;}
}void PAM::CalculateSRF()
{arma::mat tempMean(dimension, parallelSampleNumber, arma::fill::zeros);arma::mat tempVariance(dimension, parallelSampleNumber, arma::fill::zeros);for (int i = 1; i < iterationNumber + 1; ++i){for (int j = 0; j < parallelSampleNumber; ++j){tempMean.col(j) = parallelSampleMean.slice(j).col(i);tempVariance.col(j) = parallelSampleVariance.slice(j).col(i);}arma::vec B = arma::var(tempMean, 0, 1);  //对矩阵每一行求方差, N-1arma::vec W = arma::mean(tempVariance, 1); //对矩阵每一行求均值scaleReductionFactor.col(i - 1) = arma::sqrt(static_cast<double>((i - 1)) / static_cast<double>(i)+ ((parallelSampleNumber + 1) / static_cast<double>(parallelSampleNumber))* B / W);}}void PAM::PlotSRF()
{for (int i = 0; i < dimension; ++i){arma::rowvec paramSRFVec = scaleReductionFactor.row(i);  //参数均值std::vector<double> SRFVec //将arma::vec转化为std::vector<double>类型= arma::conv_to< std::vector<double> >::from(paramSRFVec);plt::plot(SRFVec);plt::title("Change process diagram of Scale Reduction Factor of parameter " + paramName.at(i));plt::xlabel("Number of simulations");plt::ylabel("Scale Reduction Factor $\\sqrt{R}$");plt::save(paramName.at(i) + "_Scale Reduction Factor" + ".png");plt::show();}
}void PAM::PlotPostDistribution()
{int size = iterationNumber - burnInNumber + 1;arma::mat tempSampleBurnIn = arma::mat(dimension, size * parallelSampleNumber, arma::fill::zeros);//去除热身期的样本for (int i = 0; i < parallelSampleNumber; ++i){tempSampleBurnIn.cols(i * size, (i + 1) * size - 1) = parallelSampleGroupBurnIn.slice(i);}for (int i = 0; i < dimension; ++i){arma::rowvec paramChainVec = tempSampleBurnIn.row(i);std::vector<double> paramChain= arma::conv_to< std::vector<double> >::from(paramChainVec);plt::hist(paramChain);plt::title("Posterior distribution histogram for parameter " + paramName.at(i));plt::xlabel(paramName.at(i));plt::ylabel("Frequency");plt::save(paramName.at(i) + "_Posterior distribution histogram" + ".png");plt::show();}}void PAM::PlotConfidenceInterval()
{CalculateWeight();SetPercentile();ecdfPred();Plot();
}void PAM::Plot()
{std::vector<double> xAxis = { -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5 };plt::plot(xAxis, sampLowerPtile);plt::plot(xAxis, sampUpperPtile);plt::plot(xAxis, sampMedian);plt::scatter(xAxis, observedData, 30, { {"color", "red"}, {"marker", "o"} });plt::title("Median and " + std::to_string(static_cast<int>(confidenceLevel * 100)) + "% AM-MCMC prediction limits with BOX");plt::xlabel("x");plt::ylabel("y");plt::save("ConfidenceInterval.png");plt::show();
}void PAM::CalculateWeight()
{int size = iterationNumber - burnInNumber + 1;int sizeTotal = iterationNumber + 1;modelledDataArray.erase(modelledDataArray.begin());NSE.erase(NSE.begin());MSE.erase(MSE.begin());BOX.erase(BOX.begin()); //删除因为PAM构造函数生成的第一行arma::vec BOXvec(BOX);  //转成arma的矢量arma::vec BOXBurnIn = arma::vec(size * parallelSampleNumber);for (int i = 0; i < parallelSampleNumber; ++i){BOXBurnIn.rows(i * size, (i + 1) * size - 1) =BOXvec.rows(burnInNumber + i * sizeTotal, iterationNumber + i * sizeTotal);}parallelBOXBurnIn= arma::conv_to< std::vector<double> >::from(BOXBurnIn);double parallelBOXBurnInSum = std::accumulate(parallelBOXBurnIn.begin(), parallelBOXBurnIn.end(), 0.0);for (auto box : parallelBOXBurnIn){double wgt = box / parallelBOXBurnInSum;weight.emplace_back(wgt);}
}void PAM::ecdfPred()
{int size = iterationNumber - burnInNumber + 1;int sizeTotal = iterationNumber + 1;modelledDataArrayBurnIn.assign(size * parallelSampleNumber, { 0 });for (int i = 0; i < parallelSampleNumber; ++i){std::copy(modelledDataArray.begin() + burnInNumber + i * sizeTotal, modelledDataArray.begin() + (i + 1) * sizeTotal, modelledDataArrayBurnIn.begin() + i * size);}modelledDataSize = modelledDataArray[0].size();for (size_t i = 0; i < modelledDataSize; ++i){std::vector<double> modValue; //某时刻所有的模拟值，不包含热身期for (auto modData : modelledDataArrayBurnIn){modValue.emplace_back(modData.at(i));}//按模拟值的升序对权重和对应的模拟值数组进行排序sort(weight, modValue);ecdf = cumsum(weight);std::vector<double> sampPtile = CalPercentiles(modValue);sampLowerPtile.emplace_back(sampPtile.at(0));sampUpperPtile.emplace_back(sampPtile.at(1));sampMedian.emplace_back(sampPtile.at(2));}
}void PAM::sort(std::vector<double>& w, std::vector<double>& x)
{//将权重与模拟值对应上，以便于扩展排序std::vector< std::pair<double, double> > wx;//构造向量wx，存储权重和模拟值对for (size_t i = 0; i < w.size(); ++i){wx.emplace_back(w[i], x[i]);}//匿名函数定义比较规则，用于pair数据结构按照权重升序排列std::sort(wx.begin(), wx.end(), [](auto a, auto b) { return a.second < b.second; });//将排序后的wx赋值给权重数组w和函模拟值数组xfor (size_t i = 0; i < wx.size(); i++){w[i] = wx[i].first;  //取出权重x[i] = wx[i].second;   //取出模拟值}
}std::vector<double> PAM::cumsum(const std::vector<double>& w)
{std::vector<double> ecdf;for (size_t i = 0; i < w.size(); i++){double wgt = std::accumulate(w.begin(), w.begin() + i, 0.0);ecdf.emplace_back(wgt);}return ecdf;
}std::vector<double> PAM::CalPercentiles(const std::vector<double>& modValue)
{std::vector<double> sampPtile;for (size_t i = 0; i < percentile.size(); i++){double percent = percentile.at(i);//构造临时向量以寻找各分位数的位置std::vector<double> tempVector;for (double val : ecdf){double temp = fabs(val - percent);tempVector.emplace_back(temp);}auto iter = std::min_element(tempVector.begin(), tempVector.end());//确保所求范围大于等于给定百分位数范围if (percent <= 0.5){if (*iter > percent){iter -= 1;}}else{if (*iter < percent){iter += 1;}}sampPtile.emplace_back(modValue.at(iter - tempVector.begin()));}return sampPtile;
}void PAM::SetPercentile()
{double lowerBound, upperBound;lowerBound = (1 - confidenceLevel) / 2.0;upperBound = 1 - lowerBound;percentile = { lowerBound, upperBound, 0.5 };
}

main.cpp

#include <iostream>#include "PAM.h"int main()
{PAM pamAlgorithm;pamAlgorithm.Estimate();}

文档说明

本地目录：E:\Research\OptA\MCMC

Gibbs：吉布斯采样估计线性回归模型参数

AM：自适应metropolis算法C++实现

RefCode:参考代码

RefPaper：参考文献

进度

2022/3/18，复现博客 Gibbs sampling for Bayesian linear regression in Python，实现吉布斯采样估计线性回归参数

2022/4/17 看MCMC首次引入水文模型参数不确定性估计的论文 Monte Carlo assessment of parameter uncertainty in conceptual catchment models Metropolis algorithm，文中所用的采样算法为Metropolis算法

2022/4/18 从Github下载AM相关代码，从知网下载AM相关论文。

论文中关于AM描述位置：

基于AM-MCMC的地震反演方法研究_李远 P27
水文模型参数优选及不确定性分析方法研究 P39
基于MCMC方法的概念性流域水文模型参数优选及不确定性研究 P63
基于AM-MCMC算法的贝叶斯概率洪水预报模型_邢贞相 P2
Haario-2001-An-adaptive-metropolis-algorithm AM算法出处

2022/4/26 阅读了上述文献中的方法介绍，准备参考E:\Research\OptA\MCMC\RefCode\adaptive_metropolis-master\adaptive_metropolis-master代码实现AM算法

2022/4/27 创建AM文件夹，用于C++实现AM算法

2022/4/28 尝试用C++实现多元正态分布抽样

2022/4/29 编写AM程序

2022/5/4 初步完成程序编写，可以运行出结果

2022/5/5 上午，完成二元一次方程测试；下午，将似然函数从纳什效率底数NSE改为均方误差MSE和BOX，效果显著，参数分布明显收敛；或许可以研究不同似然函数的影响。

2022/5/6 补充绘图、收敛判断、参考文献

2022/5/9 并行抽样

2022/5/10 完成AM-MCMC程序编写，测试线性模型参数估计通过，撰写博客

【算法】07 AM-MCMC算法C++实现相关推荐

零起点学算法07——复杂一点的表达式计算
零起点学算法07--复杂一点的表达式计算 Description 下面你来计算一个复杂一点的计算表达式 Input 没有输入 Output 输出表达式的值,保留2位小数题目分析:根号要用sqrt来算 ...
应用随机过程张波商豪_Markov链的应用一：MCMC算法
本文是张迪同学对马尔链的应用的介绍应用一:Markov链在MCMC算法中的应用 1. MCMC概念 MCMC即马尔科夫链蒙特卡洛方法(Markov Chain Monte Carlo).该方法将Ma ...
【记录读论文时遇到的一些算法1】——MCMC sampling Gibbs sampling
MCMC sampling & Gibbs sampling 1. 什么是采样 2. 为什么要采样 3. 采样的作用 4. 常用的采样方法 5. 基本采样方法 5.1 马尔科夫链蒙特卡洛采样( ...
R语言--MCMC算法介绍以及例子
Markov Chain Monte Carlo (MCMC) 是一种用于随机生成满足给定分布的样本的算法.它通过构建一个马尔可夫链来模拟一个随机过程,从而生成样本. 首先,你需要选择一个初始状态,然 ...
程序设计-在校整理-07 基于机器学习算法的DGA域名识别（NB、XGboost、MLP初探）
[在校整理-07 基于机器学习算法的DGA域名识别(NB.XGboost.MLP初探)](注:仅供参考学习使用) 一.课题内容和要求二.理论基础 2.1 DGA域名生成算法 2.2 DGA算法原理 ...
经典算法07 快速排序
经典算法07 快速排序活动地址:CSDN21天学习挑战赛 *学习的最大理由是想摆脱平庸,早一天就多一份人生的精彩:迟一天就多一天平庸的困扰. 算法思想: 在待排序表L[1-n]中任取一个元素p ...
动态规划算法-07背包问题进阶
简介我们在本专栏之前的文章介绍了基础的01背包问题及其解题思路,本文我们将讲述其拓展题型,也就是完全背包问题和多重背包问题. 01背包问题首先,我们先来简单回顾一下经典的01背包问题,关于01背包 ...
智慧交通day02-车流量检测实现07：匈牙利算法
匈牙利算法(Hungarian Algorithm)与KM算法(Kuhn-Munkres Algorithm)是用来解决多目标跟踪中的数据关联问题,匈牙利算法与KM算法都是为了求解二分图的最大匹配问题 ...
最短路径算法——Dijkstra and Floyd算法
一. 前言: 这个古老的算法应该耳熟能详了吧,但是我自从从学校出来到现在,最短路径算法都没有实际运用过,最近在一个GIS项目中总算用到了,于是乎把教材重温了下,同时查阅了网上很多的资料 ...
聚类算法当中的K-means算法如何去做天猫淘宝的推广任务
5 人赞同了该回答 figure data-size="normal">data-size="normal"> 这个入口是全网人气新品池,我们今天所 ...

【算法】07 AM-MCMC算法C++实现

目标