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

我爱C编程
280 阅读3分钟

1.算法描述

1.1先验概率的推导

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

 

1.png

 

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

 

2.png

 

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

 

3.png

4.png

 

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

 

5.png

 

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

 

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

 

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

 

6.png

 

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

 

7.png

 

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

 

8.png

 

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

 

9.png

10.png

 

2.仿真效果预览

matlab2022a仿真结果如下：

 

11.png

12.png

 

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.19dt1(i,:,kk1) - 0.02Lt(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=(2tsigma_u*(1-do_rate/t))/(D-t)/(1-(do_rate/t)/Q);

% pf -- failure pressure

pf_rate_1  = 2tsigma_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.19dt1(i,:,kk1) - 0.02Lt(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=(2tsigma_u*(1-do_rate/t))/(D-t)/(1-(do_rate/t)/Q);

    % pf -- failure pressure

    pf_rate_1  = 2tsigma_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`

avatar
文章
阅读
粉丝