基于MCMC、差分进化马尔可夫链、集成卡尔曼滤波器、近似贝叶斯计算研究（Matlab代码实现）

 👨‍🎓个人主页：研学社的博客   

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文目录如下：🎁🎁🎁

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 Matlab代码实现

---

​

💥1 概述

function DEMCD

ifferential Evolution / Differential Evolution Markov Chain. In the case of 'DE'，this function also plot a figure for the objective function of all generations.Input fields

• func : evaluation (objective) function.
• bound : 2 X n matrix, where n is the parameter dimension. Lower (first row) and upper(second row) bounds of parameter.
• size : size of population.
• life : maximum generation of population.
• funcPrior :(optional) prior distribution function to generate first generation. Default is uninformative prior within bounds.
• type :(optional) 'DE' or 'DEMC'. The former seek the global minimum, the later maps the posterior. Default is 'DEMC'.
• greedy :(optional) logical scalar, faster convergence but less stable. Breaks detailed balance and should not be used with 'DEMC'.
• logFlag :(optional) logical scalar, use log liklihood or not. Default is 'false'.
• mutation :(optional) scalar, mutation rate.
• CR :(optional) scalar, crossover rate.
• simplex :(optional) 1 X n logical maxtrix. 'true' specifies the parameters that should be positive and sum up to 1.
• name :(optional) name of output figure. Default 'DEMCtest'.

Output fields

• chain : population at each generation.
• obj : evaluation (objective) function value at each generation.
• best : best individual at last generation.

function MCMC

Markov Chain Monte Carlo with Random-Walk Metropolis.

Input fields

One of 'loc', 'funcPrior', 'bound'(in descending priority) must be specified.

• life : length of Markov Chain.
• func : liklihood function. Calculates the liklihood from parameter.
• funcProp : (optional) proposal distribution function. Make new proposal
• according to parameter. Default is normal distribution with sigma = 1.
• loc : (optional) initial parameter.
• funcPrior : (optional) prior function to generate initial parameter.
• bound : (optional) 2 X n matrix, where n is the parameter dimension.
• Lower(first row) and upper(second row) bounds of parameter.

Output fields

• chain : m X n matrix, where m = length and n = parameter dimension.
• time : m X 1 matrix, time stamp of each step.

Warning : Using 'bound' might break detailed balance and result in gibberish posterior. Not recommended.

function EnKF

Ensemble Kalman filter. Updated forecast and analysis states using observation and ensemble Kalman filter.

Input fields

• R : measurement error (sigma).
• Fstates : forecast states. 1 X n matrix, n = number of state
• Astates : analysis states. 1 X n matrix.
• obs : observation states. Scaler.
• func : model operator, from current state to next state.

Output

The same structue with updated Fstates and Astates, as well as estimation of model error (C).

function ModAv

Multi-model averaging function. Supported methods:

• EWA : Equal Weights Averaging
• BGA : Bates-Granger Averaging
• AICA : Akaike's Information Criterion Averaging
• BICA : Bayes Information Criterion Averaging
• GRA : Granger-Ramanathan Averaging
• BMA : Bayesian Model Averaging
• MMA : Mallows Model Averaging
• MMAd : Mallows Model Averaging (lie on simplex)

Input fields

• Xcali : m1 X n model prediction matrix for calibration period. m1 = time
• length of calibration period, n = model number.
• Ycali : m1 X 1 observation matrix for calibration period.
• Xeval : m2 X n model prediction matrix for evaluation period.
• Yeval : m2 X 1 observation matrix for evaluation period.
• p : 1 X n matrix of model parameter numbers.
• method : abbreviation of method.

Output fields

• weight : model weights.
• RMSEcali : RMSE of calibration period.
• RMSEeval : RMSE of evaluation period.
• chain : (BMA, MMA and MMAd) model weights in all iterations.
• obj : (BMA, MMA and MMAd) objective function value for all iterations.
• sigma : (BMA) optimized model sigma.

Warning: Method BMA could be slow and unstable.

function ABCPMC

Approximate Bayesian Computing-Population Monte Carlo method.

Input fields

• obs : 1 X d, d = dimension of target
• bound : 2 X n matrix, where n is the parameter dimension. Lower
• (first row) and upper(second row) bounds of parameter.
• funcPrior : prior distribution function to generate parameter.
• funcSummar : function to model target from parameter.
• funcDist : function to calculate distance between modeled and observed
• target.
• epsl : 1 X i epsilon matrix, where i is the number of generations.
• size : population size

Output

The same strucuture with extra fields.

• inds : individuals (parameters) for all generations.
• summar : modeled target at all generations.
• weight : weight matrix for all generations.
• C : Covariance matrix for all geneations.
• acRate : overall acceptance rate.

📚2 运行结果

​

​ 

部分代码：

%% Set default proposal distribution function to normal distribution.
if ~isfield(bead,'funcProp')
    bead.funcProp = @ (x) normrnd(x,1);
end

%% Retrive bead.loc and bead.dms (parameter dimension)
if ~isfield(bead,'loc')
    if isfield(bead,'funcPrior')
        bead.loc = bead.funcPrior();
    else
        bead.dms = size(bead.bound,2);
        bead.loc = rand(1,bead.dms)*(bead.bound(2,:)-bead.bound(1,:))+...
            bead.bound(1,:);
    end 
end
bead.dms = numel(bead.loc);

%% Create empty matrixes
rng('shuffle');
result.chain = NaN(bead.life, bead.dms);
result.time = NaN(bead.life, 1);

%% Initialize the first bead 
bead.p = bead.func(bead.loc);
bead.gen = 1;

%% Start Markov Chian and record
tic;
result.chain(bead.gen) = bead.loc;
result.time(bead.gen) = toc;
while bead.gen<=bead.life
    bead = RWM(bead);
    result.chain(bead.gen) = bead.loc;
    result.time(bead.gen) = toc;
end
end

🎉3 参考文献

部分理论来源于网络，如有侵权请联系删除。

[1]赵国强,鞠敏,刘丹丹,邓金涛,袁文文.基于卡尔曼滤波的坡度估算及验证[J].汽车电器,2022(12):68-69+72.DOI:10.13273/j.cnki.qcdq.2022.12.027.

🌈4 Matlab代码实现