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

news2025/7/9 9:04:39

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

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

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

⛳️座右铭:行百里者,半于九十。

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

目录

💥1 概述

📚2 运行结果

🎉3 参考文献

🌈4 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代码实现

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/18321.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

ShardingSphere笔记(一): 经验和踩坑总结

ShardingSphere笔记(一): 使用经验总结 文章目录ShardingSphere笔记(一): 使用经验总结一、背景框架选择二、ShardingSphere-jdbc 只是一个帮助你路由的框架(踩坑总结)1. 它默认会认…

支持末尾携带标签的多行TextView

项目开发过程中,遇到个UI上的需求,本着不重复造轮子、敏捷开发的原则,于是乎网上找寻了一番,发现还是自己搞吧,搜不到这样的需求,先看下我们的效果。 总结有以下三点需要注意: 末尾vip部分是…

$19服务:DTCStatusMask和statusofDTC bit 定义

诊断协议那些事儿 诊断协议那些事儿专栏系列文章,当ECU产生DTC时,我们只知道有故障发生了,并不清楚该故障什么时候发生,现在是否已经恢复、发生过几次,恢复过几次等信息,基于此ISO发布的14229-1使用DTC状态…

[附源码]SSM计算机毕业设计志愿者管理系统论文2022JAVA

项目运行 环境配置: Jdk1.8 Tomcat7.0 Mysql HBuilderX(Webstorm也行) Eclispe(IntelliJ IDEA,Eclispe,MyEclispe,Sts都支持)。 项目技术: SSM mybatis Maven Vue 等等组成,B/S模式 M…

【VTK+有限元后处理】可视化结果云图

构建vtkUnstructuredGrid对象 为了读取不同格式的有限元计算结果文件,我们先写一个FEDataModel类来管理有限元的几何拓扑和属性信息。 class FEDataModel:"""有限元数据模型类"""def __init__(self):self.nodes [] # 节点几何坐标…

斐波那契数列和斐波那契数

一、什么是斐波那契数列 斐波那契数列(Fibonacci sequence),又称黄金分割数列,因数学家莱昂纳多斐波那契(Leonardo Fibonacci)以兔子繁殖为例子而引入,故又称为“兔子数列”,指的是这…

【考研复试】计算机专业考研复试英语常见问题三(个人选择/学业规划篇)

相关链接: 【考研复试】计算机专业考研复试英语常见问题一(家庭/家乡/学校篇)【考研复试】计算机专业考研复试英语常见问题二(研究方向/前沿技术/本科毕设篇)【考研复试】计算机专业考研复试英语常见问题三&#xff0…

C++---哈希

目录 1. unordered系列关联式容器 1.1 unordered_map 1.1.1 unordered_map的介绍 1.1.2 unordered_map的接口说明 1.2 unordered_set 2. 底层结构 2.1 哈希概念 2.2 哈希冲突 2.3 哈希函数 2.4 哈希冲突解决 2.4.1 闭散列 2.4.2 开散列 3. 封装unorder_map和unord…

MySQL增删改查进阶 — 表的设计

文章目录表的设计1.设计思路2.实体固定关系的套路2.1 一对一关系2.2 一对多关系2.3 多对多关系3.总结表的设计 表的设计实际上要做的工作就是明确一个程序里,需要使用几个数据库,几个表,表里都有哪些列。 1.设计思路 先明确实体再明确实体…

6. Design A Web Crawler

title: Notes of System Design No.10 — Design a Web Crawler description: ‘Design a Web Crawler’ date: 2022-05-13 18:01:58 tags: 系统设计 categories: 系统设计 00. What is Web Crawler? Q :uh just for now lets just do html pagesbut your web cr…

Explaining Deepfake Detection by Analysing Image Matching 翻译

点击查看对应的代码 摘要 本文旨在解释深度伪造检测模型在仅由二进制标签做有监督时如何学习图像的伪迹特征。为此,从图像匹配的角度提出如下三个假设。1、深度伪造检测模型表明基于视觉概念的真/假图片既不与源图片相关也不与目标图片相关而是与伪迹图片相关。2、…

全链路压测效能10倍提升的压测工具实践笔记

背景 创业型公司或创新型项目往往团队资源有限,人员能力水平有限,难以投入专业自动化压测人员; 同时部分业务(tob/toc场景)长期有中小型活动场景带来小规模流量并发,需要产研能长期保障并及时感知和解决网…

GitHub Star70K登顶,字节内部数据结构与算法笔记,限时上线

为什么学算法 不得不说,现在几乎所有的大厂,比如Google、字节、BAT,面试的时候都喜欢考算法、让人现场写代码,那你有没有真正地想过,为什么这些大公司都喜欢考算法呢? 经常有人说,程序员35岁之…

Java毕业设计MVC:基于SSM实现计算机硬件评测交流平台

作者主页:编程千纸鹤 作者简介:Java、前端、Pythone开发多年,做过高程,项目经理,架构师 主要内容:Java项目开发、毕业设计开发、面试技术整理、最新技术分享 收藏点赞不迷路 关注作者有好处 项目编号&…

(杂)网易云歌单导入到apple music

喜欢apple music的简洁,就想着把网易云的歌单捣鼓进去。 获取歌单歌曲列表:https://yyrcd.com/n2s/ 转移歌单:https://soundiiz.com/zh/,首次使用需要注册,免费版只能一次导入200首。 平台选择apple music 登录授权即可…

Linux下 man命令的使用 及 中文man手册的安装

文章目录1. man命令使用2. 安装中文man手册1. man命令使用 man命令是Linux下最核心的命令之一。而man命令也并不是英文单词“man”的意思,它是单词manual的缩写,即使用手册的意思。man命令会列出一份完整的说明。其内容包括命令语法、各选项的意义及相关…

第三章 线性模型

3.1 基本形式 给定由d个属性描述的示例x(x1; x2; x3; … ; xd)。线性模型试图学得一个通过属性的线性组合来进行预测的函数,即 3.2 线性回归 线性回归试图学得一个线性模型尽可能准确地预测实值输出标记。 对于如何确定w和b,均方误差是回归任务中最常…

这样做时间轴,让你的PPT更出彩!

文章目录**▌方法一:美化时间节点****▌方法二:利用图片中的“轴”****▌方法三:时间轴不一定需要“轴”****▌方法四:把时间轴拆成数页****▌总结**已剪辑自: https://zhuanlan.zhihu.com/p/56672211 嗨,大家好&#…

【Linux】一万七千字详解 —— 基本指令(二)

文章目录前言man 指令cp 指令mv 指令echo 指令(含输出重定向)cat 指令(含输入重定向)wc 指令more 指令less 指令head 和 tail 指令(含管道用法)date 指令cal 指令sort 指令find 和 which 和 whereis 指令alias 指令grep 指令top 指令zip 和 unzip 指令结语前言 今天的主要内容…

C语言源代码系列-管理系统之学生选修课程系统

往期文章分享点击跳转>《导航贴》- Unity手册,系统实战学习点击跳转>《导航贴》- Android手册,重温移动开发 👉关于作者 众所周知,人生是一个漫长的流程,不断克服困难,不断反思前进的过程。在这个过…