clear all
close all
%BPNN-GMM.m
%Estimating GMM parameters using BPNN 
m1=2; m2=-2; v1=0.1; v2=0.2; p1=0.3; p2=0.7; 
ORIGINAL=[m1 m2 v1 v2 p1 p2]; 
x=ones(1,1000); 
for i=1:1:length(x) 
    r=rand; 
    if(r>p1) 
        t(i)=randn*sqrt(v2)+m2;
else
t(i)=randn*sqrt(v1)+m1;
end
end
[t1,x1]=hist(t,100);
t1=t1/sum(t1); 
figure(1) 
plot(x1,t1,'r')
%Number of neurons in the output layer is given %as 6 and the corresponding activation functions are given below %[m1 m2
%exp(v1) exp(v2) exp(p1)/(exp(p1)+exp(p2)) %exp(p2)/(exp(p1)+exp(p2))] the following 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% 
TRAININPUT=x';
TRAINOUTPUT=t';
W1=randn(1,100)*0.0001; 
W2=randn(6,101)*0.0001; 
JCOL=[];
for iteration=1:1:50
    YCOLTRAIN=[];
    for n=1:1:length(TRAININPUT) 
        H=logsig(TRAININPUT(n)*W1); 
        H=[H 1];
        Y=H*W2';
        %Applying activation 
        Y(5)=(exp(Y(5)))/(exp(Y(5))+exp(Y(6)));
        Y(6)=(exp(Y(6)))/(exp(Y(5))+exp(Y(6)));
        Y(3)=exp(Y(3));
        Y(4)=exp(Y(4)); 
        %Updating weights W2 
        temp1=normpdf(TRAINOUTPUT(n),Y(1),Y(3));
        temp2=normpdf(TRAINOUTPUT(n),Y(2),Y(4));
        gamma1=Y(5)*temp1/(Y(5)*temp1+Y(6)*temp2);
        gamma2=Y(6)*temp2/(Y(5)*temp1+Y(6)*temp2);
        GRADIENT=[gamma1*((Y(1)-TRAINOUTPUT(n))/(Y(3)^2)) ...
            gamma2*((Y(2)-TRAINOUTPUT(n))/(Y(4)^2))... 
            -gamma1*(((TRAINOUTPUT(n)-Y(1))^2)/(Y(3)^3)-(1/Y(3)))...
            -gamma2*(((TRAINOUTPUT(n)-Y(2))^2)/(Y(4)^3)-(1/Y(4)))... 
            (Y(5)-gamma1) (Y(6)-gamma2)]; 
        W2(1,:)=W2(1,:)-0.0001*GRADIENT(1).*H; 
        W2(2,:)=W2(2,:)-0.0001*GRADIENT(2).*H;
        W2(3,:)=W2(3,:)-0.0001*GRADIENT(3).*H;
        W2(4,:)=W2(4,:)-0.0001*GRADIENT(4).*H; 
        W2(5,:)=W2(5,:)-0.0001*GRADIENT(5).*H; 
        W2(6,:)=W2(6,:)-0.0001*GRADIENT(6).*H;
        %Updating weights W1 
        G1=GRADIENT(1).*W2(1,:).*H.*(1-H);
        G1=G1(1:1:100); 
        G2=GRADIENT(2).*W2(2,:).*H.*(1-H);
        G2=G2(1:1:100);
        G3=GRADIENT(3).*W2(3,:).*H.*(1-H);
        G3=G3(1:1:100); 
        G4=GRADIENT(4).*W2(4,:).*H.*(1-H); 
        G4=G4(1:1:100); 
        G5=GRADIENT(5).*W2(5,:).*H.*(1-H);
        G5=G5(1:1:100);
        G6=GRADIENT(6).*W2(6,:).*H.*(1-H);
        G6=G6(1:1:100);
        W1=W1-0.0001*TRAININPUT(n,:)'*(G1+G2+G3+G4+G5+G6); 
    end
    s=0; 
    for n=1:1:length(TRAININPUT)
        temp1=normpdf(TRAINOUTPUT(n),Y(1),Y(3));
        temp2=normpdf(TRAINOUTPUT(n),Y(2),Y(4));
        s=s-log(Y(5)*temp1+Y(6)*temp2);
    end
    JCOL=[JCOL s];
end
figure(2)
plot(JCOL)
title('Convergence of the Regression based BPNN')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
H=logsig(W1);
H=[H 1];
Y=H*W2';
Y(5)=exp(Y(5))/(exp(Y(5))+exp(Y(6)));
Y(6)=exp(Y(6))/(exp(Y(5))+exp(Y(6)));
Y(3)=exp(Y(3)); Y(4)=exp(Y(4));
ORIGINAL=[m1 m2 v1 v2 p1 p2];
ESTIMATED=[Y(1) Y(2) Y(3)^2 Y(4)^2 Y(5) Y(6)]; 
[ORIGINAL;ESTIMATED]
figure 
DATA=[];
for i=1:1:1000
    r=rand;
    if(r>Y(5))
        DATA=[DATA randn*Y(4)+Y(2)]; 
    else
        DATA=[DATA randn*Y(3)+Y(1)];
    end
end
[t,x]=hist(DATA,100);
t=t/sum(t);
figure(3)
plot(x,t,'k')
