数据驱动的近似贝叶斯计算框架用于流行病模型的参数估计和不确定性量化matlab复现

本文提出了一种结合动态状态识别和信息先验分布学习的贝叶斯计算框架，用于流行病模型参数估计，通过实际的COVID-19数据展示了其有效性。该方法优化了ABC算法在高维参数空间的应用，为实时流行病建模提供有力工具。

✅作者简介：热爱科研的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