✅作者简介:热爱科研的Matlab仿真开发者,修心和技术同步精进,代码获取、论文复现及科研仿真合作可私信。
🍎个人主页:Matlab科研工作室
🍊个人信条:格物致知。
更多Matlab完整代码及仿真定制内容点击👇
🔥 内容介绍
摘要
本文提出了一种用于流行病模型参数估计和不确定性量化的数据驱动的近似贝叶斯计算框架,其中包含两项新颖之处:(i) 通过使用与观测数据兼容的合理动态状态来识别初始条件;(ii) 通过交叉熵方法学习模型参数的信息先验分布。借助巴西里约热内卢市 COVID-19 流行病的实际数据,说明了新方法的有效性,采用了基于常微分方程的模型,该模型具有广义 SEIR 机制结构,包括时间相关的传播率、无症状感染者和住院人数。提出了一个具有两个成本项(住院人数和死亡人数)的最小化问题,并确定了 12 个参数。校准后的模型提供了对可用数据的持续描述,能够对未来几周的预测进行外推,这使得所提出的方法非常适用于实时流行病建模。
简介
流行病模型在预测疾病传播、评估干预措施的有效性和规划公共卫生应对措施方面发挥着至关重要的作用。然而,流行病模型通常是复杂的,具有许多未知参数。这些参数通常难以估计,因为它们可能与观测数据之间存在复杂的关系。
近似贝叶斯计算 (ABC) 是一种用于参数估计的强大方法,它可以处理复杂的模型和缺乏共轭先验分布的情况。然而,传统的 ABC 方法在计算上可能很昂贵,并且可能难以针对高维参数空间进行优化。
方法
本文提出的框架结合了两种新颖之处来克服这些挑战:(i) 通过使用与观测数据兼容的合理动态状态来识别初始条件;(ii) 通过交叉熵方法学习模型参数的信息先验分布。
初始条件识别
在流行病建模中,初始条件(例如感染者的数量)通常未知。本文提出了一种方法,通过使用与观测数据兼容的合理动态状态来识别初始条件。该方法基于以下假设:在流行病早期阶段,感染者的数量将遵循指数增长或衰减模式。
信息先验分布学习
模型参数的先验分布对于 ABC 方法的性能至关重要。本文提出了一种方法,通过交叉熵方法学习模型参数的信息先验分布。该方法基于以下思想:信息先验分布应使模型预测与观测数据尽可能接近。
正在上传…重新上传取消
📣 部分代码
% -----------------------------------------------------------------
% ABC.m
% -----------------------------------------------------------------
% -----------------------------------------------------------------
% This routine employs the ABC algorithm to solve a Bayesian
% inference problem, which seeks to find the set of parameters
% that better promotes the agreement between model predictions
% and a given set of observations.
%
% Input:
% fun - model function
% data - (Ndata x NQoI) reference data
% lb - (Nvars x 1) lower bound
% ub - (Nvars x 1) upper bound
% mu - (Nvars x 1) mean
% sigma - (Nvars x 1) standard deviation
% ABCobj - ABC object struct
%
% Output:
% X_opt - (Nvars x 1) optimal point
% F_opt - optimal value
% ABCobj - ABC object struct
%
% Reference:
% T. Toni, D. Welch, N. Strelkowa, A. Ipsen, M. Stumpf,
% Approximate Bayesian computation scheme for parameter
% inference and model selection in dynamical systems
% J. R. Soc. Interface, 6:187-202, 2007
% https://doi.org/10.1098%2Frsif.2008.0172
% -----------------------------------------------------------------
% -----------------------------------------------------------------
function [x_best,y_best,ABCobj] = ABC(fun,data,lb,ub,mu,sigma,ABCobj)
% check number of arguments
if nargin < 6
error('Too few inputs.')
elseif nargin > 7
error('Too many inputs.')
end
% check if lb or ub are empty
if isempty(lb) || isempty(ub)
error('lb and ub must be non-empty')
end
% convert lb to a column vector (if necessary)
[s1,s2] = size(lb);
if isrow(lb)
lb = lb(:);
elseif s2 ~= 1
error('lb must be a column vector')
end
% convert ub to a column vector (if necessary)
[s1,s2] = size(ub);
if isrow(ub)
ub = ub(:);
elseif s2 ~= 1
error('ub must be a column vector')
end
% check for consistency in lb and ub
if size(lb) ~= size(ub)
error('lb and ub must have the same dimensions')
end
if any(isnan(lb)) || any(~isfinite(lb))
error('lb must have finite values')
end
if any(isnan(ub)) || any(~isfinite(ub))
error('ub must have finite values')
end
%x_samples = lognrnd(mu_log,sigma_log,Nvars,Ns);
for ii=1:Nvars
x_samples(ii,:) = lognrnd(mu_log(ii),sigma_log(ii),1,Ns);
end
elseif strcmp(prior,'Gamma')
% shape and scale parameters
shape_gam = (mu./sigma).^2;
scale_gam = (sigma.^2)./mu;
% gamma
%x_samples = gamrnd(shape_gam,scale_gam,Nvars,Ns);
for ii=1:Nvars
x_samples(ii,:) = gamrnd(shape_gam(ii),scale_gam(ii),Ns,1);
end
else
% uniform
%x_samples = unifrnd(lb,ub,Nvars,Ns);
for ii=1:Nvars
x_samples(ii,:) = unifrnd(lb(ii,1),ub(ii,1),Ns,1);
end
end
% squared norm of data vector
norm_data_pow2 = data'*data;
% initiate progress bar
textprogressbar(' ');
% loop for ABC calculation
for n=1:Ns
% define model parameters vector
x = x_samples(:,n);
% compute the cost function
y = fun(x);
delta_y = data - y(1:Ndata,:);
J = sum(weigths.*diag(((delta_y')*delta_y)./norm_data_pow2));
% accept sample if the cost function is lower than a tolerance
if abs(J) < tol
% update acceptance counter
accept_counter = accept_counter + 1;
% save model parameters accepted values
ABCobj.x_accept(:,accept_counter) = x;
% save the model prediction
ABCobj.y_accept(:,accept_counter) = reshape(y,[Nfun1*Nqoi,1]);
% check if J is the minimum
if J < J_min
% save the best parameters value
x_best = x;
% save the best model prediction
y_best = y;
% update cost function minimum value
J_min = J;
end
else
% update rejection counter
reject_counter = reject_counter + 1;
% save model parameters rejected values
ABCobj.x_reject(:,reject_counter) = x;
% save the rejected model prediction
ABCobj.y_reject(:,reject_counter) = reshape(y,[Nfun1*Nqoi,1]);
end
% update progress bar
textprogressbar((n/Ns)*100);
end
% space after the progress bar
textprogressbar(' ');
fprintf('\n \n');
% update dimensions of accepted samples matrices
if accept_counter > 0
ABCobj.x_accept = ABCobj.x_accept(:,1:accept_counter);
ABCobj.y_accept = ABCobj.y_accept(:,1:accept_counter);
end
% update dimensions of rejected samples matrices
if reject_counter > 0
ABCobj.x_reject = ABCobj.x_reject(:,1:reject_counter);
ABCobj.y_reject = ABCobj.y_reject(:,1:reject_counter);
end
% number of samples
ABCobj.Ns = Ns;
% acceptance counter
ABCobj.accept_counter = accept_counter;
% rejection counter
ABCobj.reject_counter = reject_counter;
% acceptance rate
ABCobj.accept_rate = accept_counter/Ns;
% acceptance rate
ABCobj.reject_rate = reject_counter/Ns;
% error tolerance
ABCobj.tol = tol;
% minimum value for ABC error function
ABCobj.J_min = J_min;
end
% -----------------------------------------------------------------
% -----------------------------------------------------------------
% trandn
% -----------------------------------------------------------------
% This function is an efficient generator of a random vector of
% dimension length(l)=length(u) from the standard multivariate
% normal distribution, truncated over the region [l,u]. Infinite
% values for bounds 'u' and 'l' are accepted.
%
% Remark:
% If you wish to simulate a random variable 'Z' from the
% non-standard Gaussian N(m,s^2) conditional on l<Z<u, then
% first simulate X = trandn((l-m)/s,(u-m)/s) and set Z=m+s*X.
%
% Input:
% lb - (Nvars x 1) lower bound
% ub - (Nvars x 1) upper bound
%
% Output:
% x - (Nvars x 1) random vector with multiv. distribution N(0,1)
%
% References:
% Botev, Z. I. (2016). "The normal law under linear restrictions:
% simulation and estimation via minimax tilting". Journal of the
% Royal Statistical Society: Series B (Statistical Methodology).
% https://doi.org/10.1111/rssb.12162
%
% MATLAB Central File Exchange:
% Z. Botev, Truncated Normal Generator
% shorturl.at/hntuB
% -----------------------------------------------------------------
function x = trandn(l,u)
l = l(:); u = u(:); % make 'l' and 'u' column vectors
if length(l)~=length(u)
error('Truncation limits have to be vectors of the same length')
end
x = NaN(size(l));
a = .66; % treshold for switching between methods
% threshold can be tuned for maximum speed for each Matlab version
% three cases to consider:
% case 1: a < l < u
I = l > a;
if any(I)
tl = l(I); tu = u(I); x(I) = ntail(tl,tu);
end
% case 2: l < u < -a
J = u < -a;
if any(J)
t
% -----------------------------------------------------------------
⛳️ 运行结果
结果
所提出的框架在巴西里约热内卢市 COVID-19 流行病的实际数据上进行了评估。该模型能够很好地拟合观测数据,并且能够对未来几周的预测进行外推。
结论
本文提出的框架提供了一种用于流行病模型参数估计和不确定性量化的强大而有效的方法。该框架结合了两种新颖之处来克服传统 ABC 方法的挑战,并且在实际数据上显示出良好的性能。
🔗 参考文献
-
A. Cunha Jr, D. A. W. Barton, and T. G. Ritto, Uncertainty quantification in mechanistic epidemic models via cross-entropy approximate Bayesian computation, Nonlinear Dynamics, vol. 111, pp. 9649–9679, 2023 https://doi.org/10.1007/s11071-023-08327-8
Preprint available at:https://arxiv.org/abs/2207.12111
🎈 部分理论引用网络文献,若有侵权联系博主删除
🎁 关注我领取海量matlab电子书和数学建模资料
👇 私信完整代码和数据获取及论文数模仿真定制
1 各类智能优化算法改进及应用
生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化、背包问题、 风电场布局、时隙分配优化、 最佳分布式发电单元分配、多阶段管道维修、 工厂-中心-需求点三级选址问题、 应急生活物质配送中心选址、 基站选址、 道路灯柱布置、 枢纽节点部署、 输电线路台风监测装置、 集装箱船配载优化、 机组优化、 投资优化组合、云服务器组合优化、 天线线性阵列分布优化
2 机器学习和深度学习方面
2.1 bp时序、回归预测和分类
2.2 ENS声神经网络时序、回归预测和分类
2.3 SVM/CNN-SVM/LSSVM/RVM支持向量机系列时序、回归预测和分类
2.4 CNN/TCN卷积神经网络系列时序、回归预测和分类
2.5 ELM/KELM/RELM/DELM极限学习机系列时序、回归预测和分类
2.6 GRU/Bi-GRU/CNN-GRU/CNN-BiGRU门控神经网络时序、回归预测和分类
2.7 ELMAN递归神经网络时序、回归\预测和分类
2.8 LSTM/BiLSTM/CNN-LSTM/CNN-BiLSTM/长短记忆神经网络系列时序、回归预测和分类
2.9 RBF径向基神经网络时序、回归预测和分类