求解Wiener过程中随机参数分布的EM算法

嗯哼嘿嘿嘿

已于 2022-02-02 00:16:56 修改

阅读量2.9k

点赞数 9

文章标签：算法 matlab em

于 2022-02-02 00:15:05 首次发布

本文链接：https://blog.csdn.net/m0_62095064/article/details/122766994

版权

一. Wiener过程中的随机参数分布

在利用贝叶斯估方法对产品退化过程进行建模时，若以Wiener过程为模型，则其中的随机参数服从如下分布

$w=\sigma^{^{-2}}\sim Ga(a,b)$

$\mu /\omega \sim N(c,d/\omega)$

其中a,b,c,d为待求超参数。

由极大似然估计可得下式，用于M步

$\psi(\hat{a})-\ln \hat{a}=\frac{1}{M} \sum_{j=1}^{M} \ln \omega_{j}+\ln M-\ln \sum_{j=1}^{M} \omega_{j}$

$\hat{b}=\frac{M\cdot \hat{a}}{\sum_{j=1}^{M} \omega_{j}$

$\hat{c}=\frac{\sum_{j=1}^{M} \omega_{j} \mu_{j}}{\sum_{j=1}^{M} \omega_{j}}$

$\hat{d}=\frac{1}{M} \sum_{j=1}^{M}\left(\omega_{j} \mu_{j}^{2}-2 \hat{c} \omega_{j} \mu_{j}+\hat{c}^{2} \omega_{j}\right)$

上述公式中隐含数据项的期望值公式如下，用于E步

二. Matlab代码

function [para,count]=WieEM(bb)
%该函数为Winner过程的EM算法，bb=1时为直线型，bb=2时为凹凸型。
%% 初始化
[y,t]=datain;%导入性能退化数据
th=1e-6;
h=inf;
count=0;
a=0;b=0;c=0;d=0;
M=size(t,2);
ltd=cell(1,M);
yd=cell(1,M);
ydpf=cell(1,M);
s_yd=zeros(1,M);
s_ltd=zeros(1,M);
s_yltd=zeros(1,M);
n=zeros(1,M);
e1=zeros(1,M);
e2=zeros(1,M);
e3=zeros(1,M);
e4=zeros(1,M);
%% 退化数据做差
for j=1:M
    switch bb
        case 1
            for i=1:length(t{j})-1
                ltd{j}(i)=t{j}(i+1)-t{j}(i);
                yd{j}(i)=y{j}(i+1)-y{j}(i);
                ydpf{j}(i)=y{j}(i+1)^2-y{j}(i)^2;
            end
        case 2
            for i=1:length(t{j})-1
                ltd{j}(i)=t{j}(i+1)^x(3)-t{j}(i)^x(3);
                yd{j}(i)=y{j}(i+1)-y{j}(i);
                ydpf{j}(i)=y{j}(i+1)^2-y{j}(i)^2;
            end
    end
    s_yd(j)=sum(yd{j});
    s_ltd(j)=sum(ltd{j});
    s_yltd(j)=sum(ydpf{j}./2./ltd{j});
    n(j)=length(t{j});
end
%% EM迭代
while(h>=th)
    count=count+1;
    A=a;
    B=b;
    C=c;
    D=d;
    %E步
    for j=1:M
        f1=A+n(j)/2;
        f2=B+C^2/2/D-(D*s_yd(j)+C)^2/2/(D^2*s_ltd(j)+D)+s_yltd(j);
        f3=(D*s_yd(j)+C)/(D*s_ltd(j)+1);
        f4=D/(D*s_ltd(j)+1);
        e1(j)=f1/f2;
        e2(j)=psi(f1)-log(f2);
        e3(j)=e1(j)*f3;
        e4(j)=e1(j)*f3^2+f4;
    end
    E1=sum(e1);
    E2=sum(e2);
    E3=sum(e3);
    E4=sum(e4);
    %M步
    fun1=@(x)psi(x)-log(x)-E2/M-log(M)+log(E1);
    a=fsolve(fun1,1e-5);
    b=M*a/E1;
    c=E3/E1;
    d=(E4-2*c*E3+c^2*E1)/M;
    h=max([abs(A-a)/a,abs(B-b)/b,abs(C-c)/c,abs(D-d)/d]);
end
para=[a,b,c,d];
end

参考文献

[1]徐廷学,王浩伟,张鑫.EM算法在Wiener过程随机参数的超参数值估计中的应用[J].系统工程与电子技术,2015,37(03):707-712.