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