基于贝叶斯推理估计稳态 (ST) 和非稳态 (NS) LPIII 模型分布拟合到峰值放电(Matlab代码实现)

简介: 基于贝叶斯推理估计稳态 (ST) 和非稳态 (NS) LPIII 模型分布拟合到峰值放电(Matlab代码实现)

💥1 概述

在 NS 模型中,LPIII 分布的均值随时间变化。返回周期不是在真正的非平稳意义上计算的,而是针对平均值的固定值计算的。换句话说,假设分布在时间上保持固定,则计算返回周期。默认情况下,在对 NS LPIII 模型参数的估计值应用 bayes_LPIII时,将计算并显示更新的 ST 返回周期。更新的 ST 返回周期是通过计算与拟合期结束时的 NS 均值或 t = t(end) 相关的回报期获得的,假设分布在记录结束后的时间保持固定。


📚2 运行结果

部分代码:

cf  = 1;                        %multiplication factor to convert input Q to ft^3/s 
Mj  = 2;                        %1 for ST LPIII, 2 for NS LPIII with linear trend in mu
y_r = 0;                        %Regional estimate of gamma (skew coefficient)
SD_yr = 0.55;                   %Standard deviation of the regional estimate
%Prior distributions (input MATLAB abbreviation of distribution name used in 
%function call, i.e 'norm' for normal distribution as in 'normpdf')
marg_prior{1,1} = 'norm'; 
marg_prior{1,2} = 'unif'; 
marg_prior{1,3} = 'unif'; 
marg_prior{1,4} = 'unif';
%Hyper-parameters of prior distributions (input in order of use with 
%function call, i.e [mu, sigma] for normpdf(mu,sigma))
marg_par(:,1) = [y_r, SD_yr]';  %mean and std of informative prior on gamma 
marg_par(:,2) = [0, 6]';        %lower and upper bound of uniform prior on scale
marg_par(:,3) = [-10, 10]';     %lower and upper bound of uniform prior on location
marg_par(:,4) = [-0.15, 0.15]'; %lower and upper bound of uniform prior on trend 
%DREAM_(ZS) Variables
if Mj == 1; d = 3; end          %define number of parameters based on model
if Mj == 2; d = 4; end 
N = 3;                          %number of Markov chains 
T = 8000;                       %number of generations
%create function to initialize from prior
prior_draw = @(r,d)prior_rnd(r,d,marg_prior,marg_par); 
%create function to compute prior density 
prior_density = @(params)prior_pdf(params,d,marg_prior,marg_par);
%create function to compute unnormalized posterior density 
post_density = @(params)post_pdf(params,data,cf,Mj,prior_density);
%call the DREAM_ZS algorithm 
%Markov chains | post. density | archive of past states
[x,              p_x,          Z] = dream_zs(prior_draw,post_density,N,T,d,marg_prior,marg_par); 
%% Post Processing and Figures
%options:
%Which mu_t for calculating return level vs. return period? (don't change
%for estimates corresponding ti distribution at end of record, or updated ST distribution)
t = data(:,1) - min(data(:,1));                              %time (in years from the start of the fitting period)
idx_mu_n = size(t,1);                                        %calculates and plots RL vs RP for mu_t associated with t(idx_mu_n) 
                                                             %(idx_mu_n = size(t,1) for uST distribution) 
%Which return level for denisty plot? 
sRP = 100;                                                   %plots density of return level estimates for this return period
%Which return periods for output table?                      %outputs table of return level estimates for these return periods
RP_out =[200; 100; 50; 25; 10; 5; 2]; 
%end options 
%apply burn in (use only half of each chain) and rearrange chains to one sample 
x1 = x(round(T/2)+1:end,:,:);                                %burn in    
p_x1 = p_x(round(T/2)+1:end,:,:); 
post_sample = reshape(permute(x1,[1 3 2]),size(x1,1)*N,d,1); %columns are marginal posterior samples                                                           
sample_density = reshape(p_x1,size(p_x1,1)*N,1);             %corresponding unnormalized density 
%find MAP estimate of theta 
idx_theta_MAP = max(find(sample_density == max(sample_density))); 
theta_MAP = post_sample(idx_theta_MAP,:);                    %most likely parameter estimate 
%Compute mu as a function of time and credible intervals  
if Mj == 1; mu_t = repmat(post_sample(:,3),1,length(t));end  %ST model, mu is constant
if Mj == 2; mu_t = repmat(post_sample(:,3),1,length(t)) + post_sample(:,4)*t';end %NS mu = f(t)
MAP_mu_t = mu_t(idx_theta_MAP,:);                            %most likely estimate of the location parameter
low_mu_t = prctile(mu_t,2.5,1);                              %2.5 percentile of location parameter
high_mu_t = prctile(mu_t,97.5,1);                            %97.5 percentile of location parameter
%compute quantiles of the LPIII distribution 
p = 0.01:0.005:0.995;                                        %1st - 99.5th quantile (1 - 200 year RP)
a=1;
RLs = nan(size(post_sample,1),size(p,2));
for i = 1:size(post_sample,1);                               %compute return levels for each posterior sample
    RLs(i,:) = lp3inv(p,post_sample(i,1),post_sample(i,2),mu_t(i,idx_mu_n)); 
    a = a+1;
    if a == round(size(post_sample,1)/10) || i == size(post_sample,1);
    clc
    disp(['Calculating Return Levels ' num2str(round(i/size(post_sample,1)*100)) '% complete'])
    a = 1; 
    end
end
MAP_RL = RLs(idx_theta_MAP,:);                               %Return levels associated with most likely parameter estimate
low_RL = prctile(RLs,2.5,1);                                 %2.5 percentile of return level estimates
high_RL = prctile(RLs,97.5,1);                               %97.5 percentile of return level estimates


🎉3 参考文献

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

[1]Adam Luke (2022). Non-Stationary Flood Frequency Analysis .

🌈4 Matlab代码实现

相关文章
|
2天前
|
传感器 算法
ANC主动降噪理论及Matlab代码实现
ANC主动降噪理论及Matlab代码实现
车辆行驶控制运动学模型的matlab建模与仿真,仿真输出车辆动态行驶过程
该课题在MATLAB2022a中建立了车辆行驶控制运动学模型并进行仿真,展示车辆动态行驶过程。系统仿真结果包含四张图像,显示了车辆在不同时间点的位置和轨迹。核心程序定义了车辆参数和初始条件,使用ode45求解器模拟车辆运动。车辆运动学模型基于几何学,研究车辆空间位姿、速度随时间变化,假设车辆在平面运动且轮胎无滑动。运动学方程描述位置、速度和加速度关系,模型预测控制用于优化轨迹跟踪,考虑道路曲率影响,提升弯道跟踪性能。
|
30天前
|
机器学习/深度学习 算法
【MATLAB】基于VMD-SSA-LSTM的回归预测模型
【MATLAB】基于VMD-SSA-LSTM的回归预测模型
41 4
|
5天前
|
资源调度 SoC
基于UKF无迹卡尔曼滤波的电池Soc估计matlab仿真
**摘要:** 使用MATLAB2022a,基于UKF的电池SOC估计仿真比较真实值,展示非线性滤波在电动车电池管理中的效用。电池电气模型描述电压、电流与SoC的非线性关系,UKF利用无迹变换处理非线性,通过预测和更新步骤实时估计SoC,优化状态估计。尽管UKF有效,但依赖准确模型参数。
|
1月前
|
机器学习/深度学习 算法
【MATLAB】基于EMD-PCA-LSTM的回归预测模型
【MATLAB】基于EMD-PCA-LSTM的回归预测模型
35 0
【MATLAB】基于EMD-PCA-LSTM的回归预测模型
|
1月前
|
机器学习/深度学习 算法 计算机视觉
基于yolov2深度学习网络模型的鱼眼镜头中人员检测算法matlab仿真
该内容是一个关于基于YOLOv2的鱼眼镜头人员检测算法的介绍。展示了算法运行的三张效果图,使用的是matlab2022a软件。YOLOv2模型结合鱼眼镜头畸变校正技术,对鱼眼图像中的人员进行准确检测。算法流程包括图像预处理、网络前向传播、边界框预测与分类及后处理。核心程序段加载预训练的YOLOv2检测器,遍历并处理图像,检测到的目标用矩形标注显示。
|
1月前
|
算法 计算机视觉 异构计算
基于肤色模型的人脸识别FPGA实现,包含tb测试文件和MATLAB辅助验证
这是一个关于肤色检测算法的摘要:使用MATLAB 2022a和Vivado 2019.2进行测试和仿真,涉及图像预处理、RGB到YCbCr转换、肤色模型(基于阈值或概率)以及人脸检测。核心程序展示了如何读取图像数据并输入到FPGA处理,通过`tops`模块进行中值滤波、颜色空间转换及人脸检测,最终结果输出到"face.txt"。
|
1月前
|
数据安全/隐私保护
耐震时程曲线,matlab代码,自定义反应谱与地震波,优化源代码,地震波耐震时程曲线
地震波格式转换、时程转换、峰值调整、规范反应谱、计算反应谱、计算持时、生成人工波、时频域转换、数据滤波、基线校正、Arias截波、傅里叶变换、耐震时程曲线、脉冲波合成与提取、三联反应谱、地震动参数、延性反应谱、地震波缩尺、功率谱密度
基于混合整数规划的微网储能电池容量规划(matlab代码)
基于混合整数规划的微网储能电池容量规划(matlab代码)
|
1月前
|
算法 调度
含多微网租赁共享储能的配电网博弈优化调度(含matlab代码)
含多微网租赁共享储能的配电网博弈优化调度(含matlab代码)