4
$\begingroup$

I used EM algorithm to estimate the parameters of the following time invariant process. $$y=Ct+e_y$$ $$x=Pt+e_x$$

where $y \in R$, $x \in R^d$, $t\sim N(0,I^{p \times p})$, $e_y \sim N(0,\Sigma_y)$, and $e_x \sim N(0,\Sigma_x)$. Both $\Sigma_x$ and $\Sigma_y$ are diagonal covariance matrices. The set of parameters to be estimated include $\theta=(C,P,\Sigma_x, \Sigma_y)$. I first derived the joint probability density function $p(x,y,t)$ as following. $$p(x,y,t)=N(Ct,\Sigma_y)N(Pt,\Sigma_x)N(0,I)\\=x^T\Sigma_x^{-1}x+y^T\Sigma_y^{-1}y-2t^T(C^T\Sigma_y^{-1}y+P^T\Sigma_x^{-1}x)+t^T(C^T\Sigma_y^{-1}C+P^T\Sigma_x^{-1}P+I)t$$ $$$$ Knowing that $p(t|x,y)\propto p(x,y,t)$, $$Var[t|x,y]=(C^T\Sigma_y^{-1}C+P^T\Sigma_x^{-1}P+I)^{-1}$$ $$E[t|x,y]=(C^T\Sigma_y^{-1}C+P^T\Sigma_x^{-1}P+I)^{-1}(C^T\Sigma_y^{-1}y+P^T\Sigma_x^{-1}x)$$ $$Var[t|x,y]=E[(t-E[t|x,y])(t-E[t|x,y])^T]=E[tt^T|x,y]-E[t|x,y]E[t^T|x,y]$$ $$E[tt^T|x,y]=Var[t|x,y]+E[t|x,y]E[t^T|x,y]$$ $$$$ Now I have the required posterior distributions for the M-step, including $E[t] $ and $E[tt^T]$ (for conciseness $|x,y$ is omitted). In the M-step, considering the distributions for $y$ and $x$ are independent of each other, $P$, $C$, $\Sigma_x$ and $\Sigma_y$ can be maximized independent of each other as well. The complete log-likelihood function of the EM algorithm is derived as followings. $$L(X=[x_n]_{n=1}^N;\theta)=\dfrac{N}{2}\ln|\Sigma_x^{-1}|-\dfrac{1}{2}trace[\Sigma_x^{-1} \sum_{n=1}^N(x_nx_n^T-x_nE[t_n^T]P^T-PE[t_n]x_n^T+PE[t_nt_n^T]P^T)]$$

$$L(Y=[y_n]_{n=1}^N;\theta)=\dfrac{N}{2}\ln|\Sigma_y^{-1}|-\dfrac{1}{2}trace[\Sigma_y^{-1} \sum_{n=1}^N(y_ny_n^T-y_nE[t_n^T]C^T-CE[t_n]y_n^T+CE[t_nt_n^T]C^T)]$$

Subsequently, taking the derivatives of the above equations with respect to each parameter and setting these derivatives to zero yield (derivation similar to equations 13 to 16 of http://mlg.eng.cam.ac.uk/zoubin/course04/tr-96-2.pdf): $$C_{new}=(\sum_{n=1}^Ny_nE[t_n^T])(\sum_{n=1}^NE[t_nt_n^T])^{-1}$$ $$P_{new}=(\sum_{n=1}^Nx_nE[t_n^T])(\sum_{n=1}^NE[t_nt_n^T])^{-1}$$ $$\Sigma_x=\dfrac{1}{N}\sum_{n=1}^N(x_nx_n^T-P_{new}E[t_n]x_n^T)$$ $$\Sigma_y=\dfrac{1}{N}\sum_{n=1}^N(y_ny_n^T-C_{new}E[t_n]y_n^T)$$ This completes one step of EM. The joint likelihood of the data is computed as following, since $p(x)$ and $p(y)$ are independent of each other (not sure if it is correct).This log-likelihood should increase after each step of EM. $$L(X,Y;\theta_{new})=N\ln(2\pi)^{-1/2}+N\ln(2\pi)^{-d/2}+N\ln|\Sigma_y+C_{new}C^T_{new}|^{-1/2}+N\ln|\Sigma_x+P_{new}P^T_{new}|^{-1/2}-\dfrac{1}{2}\sum_{n=1}^Nx_n^T(\Sigma_x+P_{new}P^T_{new})^{-1}x_n+y_n^T(\Sigma_y+C_{new}C^T_{new})^{-1}y_n$$

I have implemented this in matlab. The log-likelihood is decreasing after each step. I am attaching the code and the testing data with this question. Could anybody help me identify the problem?

function net=PPLS_EM(X,Y,p,iter,tol)

net=struct('type','PPLS_EM','C',[],'P',[],'Ex',[],'Ey',[],'mux',[],'muy',    [],'LL',[]);
[N,din]=size(X); %N by d input data matrix column [6 9 11 13 16 18] of the d00_te. data
[~,dout]=size(Y); %N by 1 output data matrix column 7 of the d00_te.dat

const1=(2*pi)^(-din/2);
const2=(2*pi)^(-dout/2);
lik=0;
LL=[];
tiny=exp(-700);

Vcur=zeros(p,p,N);
tcur=zeros(p,N);

YY=diag(sum(Y.*Y)');
XX=diag(sum(X.*X)');

if nargin<5
tol=0.0001;
end

if nargin<4
iter=100;
end

if nargin<3
p=2;
end


%% Initilization
C=ones(dout,p);
Ey=ones(dout,1);
Ey=diag(Ey);

P=ones(din,p);
Ex=ones(din,1);
Ex=diag(Ex);

Ex=Ex+(Ex==0)*tiny;
Ey=Ey+(Ey==0)*tiny;

Etxsum=0;
Etysum=0;
Vttsum=0;

for t=1:iter
%% E-step

oldlik=lik;
for n=1:N
    tcur(:,n)=inv(eye(p)+P'*inv(Ex)*P+C'*inv(Ey)*C)*(P'*inv(Ex)*X(n,:)'+C'*inv(Ey)*Y(n,:)');
    Vcur(:,:,n)=inv(eye(p)+P'*inv(Ex)*P+C'*inv(Ey)*C)+tcur(:,n)*tcur(:,n)';

    invx=inv(Ex+P*P');
    invy=inv(Ey+C*C');

    Etx=tcur(:,n)*X(n,:);
    Ety=tcur(:,n)*Y(n,:);

    Etxsum=Etxsum+Etx;
    Etysum=Etysum+Ety;

    Vttsum=Vttsum+Vcur(:,:,n);

    % calculate likelihood
    detiEx=sqrt(det(invx));
    detiEy=sqrt(det(invy));
    if (isreal(detiEx) && detiEx>0)
        lik=lik+log(detiEx)-0.5*sum(sum(X(n,:).*(X(n,:)*invx)))+log(detiEy)-0.5*sum(sum(Y(n,:).*(Y(n,:)*invy))); %log-likelihood of input data matrix
    else
        break;
    end;
end
lik=lik+N*log(const1)+N*log(const2);

if (t<=2)
    likbase=lik;
elseif (lik<oldlik)
    %fprintf('Oops!');
elseif ((lik-likbase)<(1 + tol)*(oldlik-likbase)||~isfinite(lik))
    break;
end;

LL=[LL lik];

%% M-step
C=Etysum'*inv(Vttsum);
P=Etxsum'*inv(Vttsum);
Ex=diag(diag(XX-(P*Etxsum)))/N;
Ey=diag(diag(YY-(C*Etysum)))/N;
end

net.C=C;
net.P=P;
net.Ex=Ex;
net.Ey=Ey;
net.LL=LL;

Data can be accessed at https://www.dropbox.com/s/i36kw4qq6al0pmn/d00_te.dat?dl=0

Results

$C_{new}$ for $p=2$ for 10 steps val(:,:,1) =

-6.3224    1.7853

val(:,:,2) =

-6.3503    1.7866

val(:,:,3) =

-2.9558    1.9637

val(:,:,4) =

-1.6677    1.0381

val(:,:,5) =

-1.6664    1.0375

val(:,:,6) =

-1.6651    1.0368

val(:,:,7) =

-1.6637    1.0362

val(:,:,8) =

-1.6623    1.0355

val(:,:,9) =

-1.6609    1.0349

val(:,:,10) =

-1.6595    1.0342

$P_{new}$ for 10 steps

val(:,:,1) =

 0.0097   -0.0013
 0.0030    0.0199
 0.1186    0.0240
-6.5193    1.8837
-5.6963    1.8107
 0.2769   -0.0450

val(:,:,2) =

 0.0098   -0.0013
 0.0030    0.0198
 0.1193    0.0237
-6.5530    1.8852
-5.6773    1.8080
 0.2731   -0.0448

val(:,:,3) =

 0.0052    0.0015
-0.0059    0.0052
 0.0591   -0.0072
-3.0491    2.0738
-2.6955    1.7138
 0.1382   -0.0148

val(:,:,4) =

 0.0042    0.0028
-0.0050    0.0045
 0.0338    0.0123
-1.7145    1.1179
-1.5561    0.8878
 0.0927    0.0243

val(:,:,5) =

 0.0042    0.0028
-0.0050    0.0045
 0.0338    0.0123
-1.7131    1.1172
-1.5549    0.8872
 0.0927    0.0244

val(:,:,6) =

 0.0042    0.0028
-0.0050    0.0045
 0.0338    0.0123
-1.7116    1.1165
-1.5537    0.8867
 0.0927    0.0244

val(:,:,7) =

 0.0042    0.0028
-0.0050    0.0045
 0.0338    0.0123
-1.7101    1.1158
-1.5526    0.8861
 0.0927    0.0244

val(:,:,8) =

 0.0042    0.0027
-0.0050    0.0045
 0.0338    0.0123
-1.7086    1.1151
-1.5514    0.8855
 0.0926    0.0244

val(:,:,9) =

 0.0043    0.0027
-0.0049    0.0045
 0.0337    0.0124
-1.7071    1.1144
-1.5501    0.8849
 0.0926    0.0244

val(:,:,10) =

 0.0043    0.0027
-0.0049    0.0045
 0.0337    0.0124
-1.7055    1.1136
-1.5489    0.8843
 0.0926    0.0244

Log likelihood for 10 steps -6176.77412636697 + 0.00000000000000i
-12350.0028284986 + 0.00000000000000i
-48330.2798829560 + 0.00000000000000i
-54546.8694822146 + 1507.96447372310i
-60722.1364253500 + 1507.96447372310i
-66897.4033684854 + 1507.96447372310i
-73072.6703116208 + 1507.96447372310i
-79247.9372547562 + 1507.96447372310i
-85423.2041978916 + 1507.96447372310i
-91598.4711410270 + 1507.96447372310i

Improve this question
$\endgroup$
3
  • $\begingroup$ Could you provide a log of both the new parameters and the new loglikelihood after each step (lets say for the first 10 steps) ? $\endgroup$
    – Nikolas Rieble
    Commented Dec 11, 2016 at 13:18
  • $\begingroup$ $Cnew$ for $p=2$ for 10 steps val(:,:,1) = -6.3224 1.7853 val(:,:,2) = -6.3503 1.7866 val(:,:,3) = -2.9558 1.9637 val(:,:,4) = -1.6677 1.0381 val(:,:,5) = -1.6664 1.0375 val(:,:,6) = -1.6651 1.0368 val(:,:,7) = -1.6637 1.0362 val(:,:,8) = -1.6623 1.0355 val(:,:,9) = -1.6609 1.0349 val(:,:,10) = -1.6595 1.0342 $\endgroup$
    – Karlsson Yu
    Commented Dec 11, 2016 at 13:59
  • $\begingroup$ $Pnew$ for 3 steps too long to be included val(:,:,1) = 0.0097 -0.0013 0.0030 0.0199 0.1186 0.0240 -6.5193 1.8837 -5.6963 1.8107 0.2769 -0.0450 val(:,:,2) = 0.0098 -0.0013 0.0030 0.0198 0.1193 0.0237 -6.5530 1.8852 -5.6773 1.8080 0.2731 -0.0448 val(:,:,3) = 0.0052 0.0015 -0.0059 0.0052 0.0591 -0.0072 -3.0491 2.0738 -2.6955 1.7138 0.1382 -0.0148 $\endgroup$
    – Karlsson Yu
    Commented Dec 11, 2016 at 14:00

1 Answer 1

Reset to default
1
$\begingroup$

Problem Solved

The decreasing Likelihood was due to a coding problem. The accumulation of $lik$ is not only in the for loop for the data samples but also in the for loop for E-M step. This leads to adding the $lik$ of current EM-step onto the $lik$ of the previous EM-step which of course keeps decreasing $lik$. I am attaching the revised code here for anyone interested. No restriction on how you might use it.

function net=PPLS_EM(X,Y,p,iter,tol)

net=struct('type','PPLS_EM','C',[],'P',[],'Ex',[],'Ey',[],'LL',[]);
[N,~]=size(X); %N by p input data matrix

LL=[];
lik=0;

YY=diag(sum(Y.*Y)');
XX=diag(sum(X.*X)');

if nargin<5
tol=0.0001;
end

if nargin<4
iter=100;
end

if nargin<3
p=2;
end


%% Initilization

C=ones(dout,p);
Ey=ones(dout,1);
Ey=diag(Ey);

P=ones(din,p);
Ex=ones(din,1);
Ex=diag(Ex);

for t=1:iter
%% E-step
oldlik=lik;
[lik,Vtt,Etxsum,Etysum]=PPLS_E(C,P,Ex,Ey,X,Y);
if (t<=2)
    likbase=lik;
elseif (lik<oldlik)
    %fprintf('Oops!');
elseif ((lik-likbase)<(1 + tol)*(oldlik-likbase)||~isfinite(lik))
    break;
end;

LL=[LL lik];

%% M-step
C=Etysum'*inv(Vtt);
P=Etxsum'*inv(Vtt);
Ex=diag(diag(XX-(P*Etxsum)))/N;
Ey=diag(diag(YY-(C*Etysum)))/N;
end

net.C=C;
net.P=P;
net.Ex=Ex;
net.Ey=Ey;
net.LL=LL;

And the code for PPLS_E

    function [lik,Vtt,Etxsum,Etysum]=PPLS_E(C,P,Ex,Ey,X,Y)

    [N,din]=size(X);
    [~,dout]=size(Y);

    const1=(2*pi)^(-din/2);
    const2=(2*pi)^(-dout/2);

    lik=0;
    p=size(P,2);

    Vcur=zeros(p,p,N);
    tcur=zeros(p,N);

    tiny=exp(-700);
    I=eye(p);

    Ex=Ex+(Ex==0)*tiny;
    Ey=Ey+(Ey==0)*tiny;

    Etxsum=zeros(p,din);
    Etysum=zeros(p,dout);
    Vtt=0;

    Sigma_XY=I+C'*inv(Ey)*C+P'*inv(Ex)*P;
    for n=1:N
        tcur(:,n)=inv(Sigma_XY)*(P'*inv(Ex)*X(n,:)'+C'*inv(Ey)*Y(n,:)');
        Vcur(:,:,n)=inv(Sigma_XY)+tcur(:,n)*tcur(:,n)';

        Etx=tcur(:,n)*X(n,:);
        Ety=tcur(:,n)*Y(n,:);

        Etxsum=Etxsum+Etx;
        Etysum=Etysum+Ety;

        Vtt=Vtt+Vcur(:,:,n);

        % calculate likelihood
        Ydiff=Y(n,:)'-C*tcur(:,n);
        Xdiff=X(n,:)'-P*tcur(:,n);
        inx=inv(Ex+P*Vcur(:,:,n)*P');
        iny=inv(Ey+C*Vcur(:,:,n)*C');
        detiEx=sqrt(det(inx));
        detiEy=sqrt(det(iny));
        if (isreal(detiEx) && detiEx>0 && isreal(detiEy) && detiEy>0)
            lik=lik+log(detiEx)-0.5*sum(Xdiff'.*(Xdiff'*inx))+log(detiEy)-0.5*sum(Ydiff'.*(Ydiff'*iny)); %log-likelihood of input data matrix
        else
            break;
        end;
    end
    lik=lik+N*log(const1)+N*log(const2);
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.