首页 > 其他分享 >基于概率论的MATLAB仿真,内容包括非共轭条件下的后验概率的推导,共轭条件下的非完备集的后验概率的推导

基于概率论的MATLAB仿真,内容包括非共轭条件下的后验概率的推导,共轭条件下的非完备集的后验概率的推导

时间:2023-03-05 13:34:08浏览次数:53  
标签:do 后验 推导 kk1 rate pf dt1 共轭

1.算法描述

1.1先验概率的推导

        根据贝叶斯概率论可知,某一事件的后验概率可以根据先验概率来获得,因此,这里首先对事件的先验概率分布进行理论的推导。假设测量的腐蚀数据服从gamma分布,其概率密度函数可以通过如下表达式表示:

 

 

 

 

       根据参考文献1和参考文献2的理论推导可知,采用反gamma分布,可以作为腐蚀数据的先验分布,即:

 

 

 

 

       公式3为公式2的自然指数形式,公式3中,x表示腐蚀数据,参数a和b分别表示反gamma分布的参数估计值。

 

 

 

 

 

 

从公式7可知,此时后验概率值则取决于最后一次测量结果.根据上述推导过程,完备集的后验概率可以通过如下公式计算得到:

 

 

 

 

       但是完备集下的后验概率所满足的公式3条件和公式4条件,在实际中往往不太可能发生,因此需要考虑非完备集下的后验概率计算方法。  

 

1.2.共轭条件下的非完备集的后验概率的推导

 

        完备集下的后验概率不太适用于实际情况,因此,对于实际情况,需要考虑非完备集下的后验概率的计算。非完备集下的后验概率是关于随机事件的条件概率,是在相关证据给定并纳入考虑之后的条件概率。后验概率和先验概率满足如下关系式:

 

 

 

 

      从公式可知,后验概率等同于先验函数和似然函数的乘积,这里先验函数根据本文公式2获得,下面主要对似然函数进行公式推导,根据参考文献5的相关推导过程可知,后验概率的基本计算公式如下:  

 

 

 

 

根据本文上述章节的介绍,参数A和B满足如下关系式:

 

 

 

 

因此,似然函数可以通过如下表达式表示:

 

 

 

 

 

 

2.仿真效果预览

matlab2022a仿真结果如下:

 

 

 

 

 

 

3.MATLAB核心程序

 

K_d        = length(dt(:,:,kk1)); %total number of d
K_l        = length(Lt(:,:,kk1)); %total number of l
 
for i = 1:K_d
    if Nn2(i) == 1
       dt1(i,:,kk1) = dt1(i,:,kk1); 
    else
       dt1(i,:,kk1) = 5.39 + 0.19*dt1(i,:,kk1) - 0.02*Lt(i,:,kk1) + 0.35*Nn2(i);
    end
end
%m->mm
dt1        = 1000*dt1;
%to obtaion a average number of do_rate and Lo_rate
do_rate    = sum(dt1(:,:,kk1))/K_d;  
Lo_rate    = sum(Lt(:,:,kk1))/K_l; 
% Q = sqrt(1+0.31*power(Lo_rate/sqrt(D/t),2)); 
% Q--length of correction factor
Q1         =(Lo_rate/sqrt(D_t))^2;
Q          = sqrt(1+0.31*Q1);
% pf_rate=(2*t*sigma_u*(1-do_rate/t))/(D-t)/(1-(do_rate/t)/Q);
% pf -- failure pressure
pf_rate_1  = 2*t*sigma_u*(1-do_rate/t);
pf_rate_2  =(D-t)*(1-do_rate/t/Q);
pf_rate    = pf_rate_1/pf_rate_2;
grid_dist  = 0.1/20; % in order to get the obvious result on the plot
x          = grid_dist:grid_dist:pf_rate*0.015;
%fit the contineous inverted gamma density to the data
par        = invgamafit(0.1); % change pf_rate from mPa to kPa, in order to get the obvious result on the plot
a          = par(1);
b          = 1/par(2);
%Examining inverted gamma distributed prior
prior     = exp(a*log(b)-gammaln(a)+(-a-1)*log(x)-b./x);
load r2.mat
prior     = post_imp_prior';
%Examination of inverted gamma post prior after perfect inspection
A         = a + dt1(K_d)/pf_rate^2;
B         = b +  Lt(K_l)/pf_rate^2;
postprior = exp(A*log(B)-gammaln(A)-(A+1)*log(x)-B./x);
%***********************************************************************************
% % %***********************************************************************************
% %定义likelyhood
% likeliprod = likelihoods(x,t,dt(:,:,kk1),Lt(:,:,kk1),Nn2);
%***********************************************************************************
%这个部分和之前的不一样了,修改后的如下所示:
%***********************************************************************************
%对prior参数进行随机化构造
m = 10;
for ijk = 1:m
    ijk
    %***********************************************************************************
    %***********************************************************************************
    %Calaulate the depth change rate and length change rate with time 
    for kk1 =1:(kk -1);
        drate1 = normrnd(drate,drateS, nsamples,1, kk1); % Measured defect depth @ time T 
        Lrate1 = normrnd(Lrate,LrateS, nsamples,1, kk1); % Measured defect length @ time T    
        if kk1 == 1
           dt(:,:,kk1) = do1(:,:,kk1) + drate1(:,:,kk1)*(delT) ; 
           dt1(:,:,kk1) = dt(:,:,kk1);
           Lt(:,:,kk1) = Lo1(:,:,kk1) + Lrate1(:,:,kk1)*(delT) ;    
        else 
           dt(:,:,kk1) = dt(:,:,kk1-1)   + drate1(:,:,kk1)*(delT);
           dt1(:,:,kk1) = dt(:,:,kk1) ;
           Lt(:,:,kk1) = Lt(:,:,kk1-1) + Lrate1(:,:,kk1)*(delT); 
        end  
    end
    K_d        = length(dt(:,:,kk1)); %total number of d
    K_l        = length(Lt(:,:,kk1)); %total number of l
    for i = 1:K_d
        if Nn2(i) == 1
           dt1(i,:,kk1) = dt1(i,:,kk1); 
        else
           dt1(i,:,kk1) = 5.39 + 0.19*dt1(i,:,kk1) - 0.02*Lt(i,:,kk1) + 0.35*Nn2(i);
        end
    end
    %m->mm
    dt1        = 1000*dt1;
    %to obtaion a average number of do_rate and Lo_rate
    do_rate    = sum(dt1(:,:,kk1))/K_d;  
    Lo_rate    = sum(Lt(:,:,kk1))/K_l; 
    % Q = sqrt(1+0.31*power(Lo_rate/sqrt(D/t),2)); 
    % Q--length of correction factor
    Q1         =(Lo_rate/sqrt(D_t))^2;
    Q          = sqrt(1+0.31*Q1);
    % pf_rate=(2*t*sigma_u*(1-do_rate/t))/(D-t)/(1-(do_rate/t)/Q);
    % pf -- failure pressure
    pf_rate_1  = 2*t*sigma_u*(1-do_rate/t);
    pf_rate_2  =(D-t)*(1-do_rate/t/Q);
    pf_rate    = pf_rate_1/pf_rate_2;
    grid_dist  = 0.1/20; % in order to get the obvious result on the plot
    x          = grid_dist:grid_dist:pf_rate*0.015;
    %fit the contineous inverted gamma density to the data
    par        = invgamafit(0.1); % change pf_rate from mPa to kPa, in order to get the obvious result on the plot
    as(1,ijk)  = par(1);
    bs(1,ijk)  = 1/par(2);
    %***********************************************************************************
    %***********************************************************************************
end

 

  

 

标签:do,后验,推导,kk1,rate,pf,dt1,共轭
From: https://www.cnblogs.com/51matlab/p/17180318.html

相关文章

  • 基于概率论的MATLAB仿真,内容包括非共轭条件下的后验概率的推导,共轭条件下的非完备集
    1.算法描述1.1先验概率的推导根据贝叶斯概率论可知,某一事件的后验概率可以根据先验概率来获得,因此,这里首先对事件的先验概率分布进行理论的推导。假设测量的腐蚀数据服从g......
  • pearson总体相关系数到样本相关系数推导过程
    相关系数是根据样本数据计算的度量两个变量之间线性关系强的统计量,若相关关系是根据总体全部数据计算的,成为总体相关系数,记为\(\rho\);若是根据样本数据计算的,则称为样本相......
  • 推导式
    """应用场景:循环次数固定的遍历:从头到尾全部访问一遍对象:列表/字典#把列表/字典每一个成员,取出来,赋值给变量for变量名(随便取)in列表/字典:取到的第一个成员......
  • 高斯判别分析GDA推导与代码实现
    高斯判别分析GDA推导与代码实现生成学习处理分类问题,我们可以使用逻辑回归、Softmax。这两种方法都属于“判别学习”,也就是给定\((x^{(i)},y^{(i)})\),我们学习\(P(y|x......
  • python各种推导式(comprehensions)
    各种推导式(comprehensions)推导式(又称解析式)是Python的一种独有特性,如果我被迫离开了它,我会非常想念。推导式是可以从一个数据序列构建另一个新的数据序列的结构体。共......
  • 圆锥曲线 / conic section の 推导 / proof
    #####椭圆/eclipse-equation/公式:$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1\/\\frac{y^2}{a^2}+\frac{x^2}{b^2}=1\quad\left(a>b>0\right)$-顶点:$(\pma,0),(0,\pmb......
  • 单层感知机推导?
    1.单层感知器推导 2.推导公式3.总结 TRANSLATEwithxEnglishArabicHebrewPolishBulgarianHindiPortugueseCata......
  • 期权定价之希腊值推导(Derivation of Greeks in European Option Pricing)
    DerivationofGreeksinEuropeanOptionPricingTheEuropeanVanillacallandputoptionswithdividendpaymentshavethevaluesinformulaswith:\[\begin{al......
  • 第一篇:Web框架推导
    纯手工撸web框架importsocketserver=socket.socket()server.bind(('127.0.0.1',8080))server.listen(5)whileTrue:conn,addr=server.accept()dat......
  • python列表推导式的结构探究
    1、列表推导式结构包含在一对方括号中,一个表达式,后面是for子句,然后是零个或多个for或if子句。2、其结果将是一个新列表,根据for和if子句的内容计算表达式。实例fromcollecti......