%MODEL SPECIFICATIONS
Times=input('Give the number of time points')
J=input('Give the number of clusters at each time point in vector form ')
T=input('Give the number of iterations T  ')
F=input('Give the number of iterations to throw away in the beginning (Burn-in) F  ')
Default=input('To use default prior specification, type 0. To specify own priors, type 1  ')

[D(1),n] = size(Y1);     
if D(1)>n
    warning('The number of observations n have to exceed the dimension K of data. Respecify the data matrix Y1')
    break
end

[D(2),n] = size(Y2);
if D(2)>n
    warning('The number of observations n have to exceed the dimension K of data. Respecify the data matrix Y2')
    break
end

if Times==3
    Y3=Y3'
    [D(3),n] = size(Y3);
    if D(3)>n
    warning('The number of observations n have to exceed the dimension K of data. Respecify the data matrix Y2')
    break
    end
end

%DEFAULT PRIOR SPECIFICATION
if Default==0 
    
    PRIOR1tau=ones(1,J(1));
    PRIOR2tau=ones(1,J(2));
    
    PRIOR1ksi=repmat(mean(Y1')',1,J(1));
    PRIOR2ksi=repmat(mean(Y2')',1,J(2));

    PRIOR1m=10*ones(1,J(1));
    PRIOR2m=10*ones(1,J(2));
    
    Covariance1=cov(Y1');
    Covariance2=cov(Y2');
    for i=1:J(1)
        PRIOR1.sig(:,i*D(1)-(D(1)-1):i*D(1))=Covariance1;
        PRIOR1psi(:,i*D(1)-(D(1)-1):i*D(1))=PRIOR1m(i)*Covariance1;
    end
    for i=1:J(2)
        PRIOR2.sig(:,i*D(2)-(D(2)-1):i*D(2))=Covariance2;
        PRIOR2psi(:,i*D(2)-(D(2)-1):i*D(2))=PRIOR2m(i)*Covariance2;
    end
    PRIOR1.psi=PRIOR1psi
    PRIOR2.psi=PRIOR2psi

    PRIOR1alfa=5*ones(1,J(1));
    PRIOR2alfa=5*ones(1,J(2));
    
    PRIOR1beta=5*ones(J(1),J(2));
   
    if Times==3
        PRIOR3tau=ones(1,J(3));
        PRIOR3ksi=repmat(mean(Y3')',1,J(3));
        PRIOR3m=10*ones(1,J(3));
        Covariance3=cov(Y3');
        for i=1:J(2)
        PRIOR3.sig(:,i*D(3)-(D(3)-1):i*D(3))=Covariance3;
        PRIOR3psi(:,i*D(3)-(D(3)-1):i*D(3))=PRIOR3m(i)*Covariance3;
        end
        PRIOR3.psi=PRIOR3psi
        PRIOR2beta=5*ones(J(2),J(3));
    end
        
else
    
% MANUAL PRIOR SPECIFICATINS
    %PRIOR SPECIFICATIONS - CLUSTER MEANS
    %Time1
    PRIOR1tau=input('According to Instruction 1, Section 3.4, specify the precision parameters for the mean variables at Time 1  ') 
    %WARNINGS
    [Ro,Co]=size(PRIOR1tau);
    while Ro~=1
        warning('This is a vector, i.e. the number of rows has to be equal to 1')
        PRIOR1tau=input('According to Instruction 1, Section 3.4, retype the precision parameters for the mean variables at Time 1  ') 
        [Ro,Co]=size(PRIOR1tau);
    end
    while Co~=J(1)
        warning('The number of columns has to equal the number of clusters at time 1')
        PRIOR1tau=input('According to Instruction 1, Section 3.4, retype the precision parameters for the mean variables at Time 1  ') 
        [Ro,Co]=size(PRIOR1tau);
    end
    
    %Time 2
    PRIOR2tau=input('According to Instruction 1, Section 3.4, specify the precision parameters for the mean variables at Time 2  ') 
    %WARNINGS
    [Ro,Co]=size(PRIOR2tau);
    while Ro~=1
        warning('This is a vector, i.e. the number of rows has to be equal to 1')
        PRIOR2tau=input('According to Instruction 1, Section 3.4, retype the precision parameters for the mean variables at Time 2  ') 
        [Ro,Co]=size(PRIOR2tau);
    end
    while Co~=J(2)
        warning('The number of columns has to equal the number of clusters at time 2')
        PRIOR2tau=input('According to Instruction 1, Section 3.4, retype the precision parameters for the mean variables at Time 2  ') 
        [Ro,Co]=size(PRIOR2tau);
    end
    
    %Time 3
    if Times==3
    PRIOR3tau=input('According to Instruction 1, Section 3.4, specify the precision parameters for the mean variables at Time 3  ') 
    %WARNINGS
    [Ro,Co]=size(PRIOR3tau);
    while Ro~=1
        warning('This is a vector, i.e. the number of rows has to be equal to 1')
        PRIOR3tau=input('According to Instruction 1, Section 3.4, retype the precision parameters for the mean variables at Time 3  ') 
        [Ro,Co]=size(PRIOR3tau);
    end
    while Co~=J(3)
        warning('The number of columns has to equal the number of clusters at time 3')
        PRIOR3tau=input('According to Instruction 1, Section 3.4, retype the precision parameters for the mean variables at Time 3  ') 
        [Ro,Co]=size(PRIOR3tau);
    end
    end

 
    %Time 1
    PRIOR1ksi=input('According to Intruction 2, Section 3.4, specify the prior of the mean values for all clusters at Time 1  ')
    %WARNINGS
    [Ro,Co]=size(PRIOR1ksi);
    while Ro~=D(1)
        warning('The number of rows has to equal the number of variables D(1)')
        PRIOR1ksi=input('According to Intruction 2, Section 3.4, retype the prior of the mean values for all clusters at Time 1  ')
        [Ro,Co]=size(PRIOR1ksi);
    end
    while Co~=J(1)
        warning('The number of columns has to equal the number of clusters at time 1')
        PRIOR1ksi=input('According to Intruction 2, Section 3.4, retype the prior of the mean values for all clusters at Time 1  ')
        [Ro,Co]=size(PRIOR1ksi);
    end
    
    %Time 2
    PRIOR2ksi=input('According to Intruction 2, Section 3.4, specify the prior of the mean values for all clusters at Time 2  ')
    %WARNINGS
    [Ro,Co]=size(PRIOR2ksi);
    while Ro~=D(2)
        warning('The number of rows has to equal the number of variables D(2)')
        PRIOR2ksi=input('According to Intruction 2, Section 3.4, retype the prior of the mean values for all clusters at Time 2  ')
        [Ro,Co]=size(PRIOR2ksi);
    end
    while Co~=J(2)
        warning('The number of columns has to equal the number of clusters at time 2')
        PRIOR2ksi=input('According to Intruction 2, Section 3.4, retype the prior of the mean values for all clusters at Time 2  ')
        [Ro,Co]=size(PRIOR2ksi);
    end

    %Time 3
    if Times==3
    PRIOR3ksi=input('According to Intruction 2, Section 3.4, specify the prior of the mean values for all clusters at Time 3  ')
    %WARNINGS
    [Ro,Co]=size(PRIOR3ksi);
    while Ro~=D(3)
        warning('The number of rows has to equal the number of variables D(3)')
        PRIOR3ksi=input('According to Intruction 2, Section 3.4, retype the prior of the mean values for all clusters at Time 3  ')
        [Ro,Co]=size(PRIOR3ksi);
    end
    while Co~=J(3)
        warning('The number of columns has to equal the number of clusters at time 3')
        PRIOR3ksi=input('According to Intruction 2, Section 3.4, retype the prior of the mean values for all clusters at Time 3  ')
        [Ro,Co]=size(PRIOR3ksi);
    end
    end

    %PRIOR SPECIFICATIONS - CLUSTER COVARIANCES
    %Time 1
    PRIOR1m=input('According to Instruction 3, Section 3.4, specify the degrees of freedom for the covariance matrices at Time 1    ')
    %WARNINGS
    [Ro,Co]=size(PRIOR1m);
    while Ro~=1
        warning('This is a vector, i.e. the number of rows has to be equal to 1')
        PRIOR1m=input('According to Instructions 3 Section 3.4, retype the degrees of freedom for the covariance matrices at Time 1   ')
        [Ro,Co]=size(PRIOR1m);
    end
    while Co~=J(1)
        warning('The number of columns has to equal the number of clusters at time 1')
        PRIOR1m=input('According to Instructions 3, Section 3.4, retype the degrees of freedom for the covariance matrices at Time 1   ')
        [Ro,Co]=size(PRIOR1m);
    end
    
    %Time 2
    PRIOR2m=input('According to Instruction 3, Section 3.4, specify the degrees of freedom for the covariance matrices at Time 2    ')
    %WARNINGS
    [Ro,Co]=size(PRIOR2m);
    while Ro~=1
        warning('This is a vector, i.e. the number of rows has to be equal to 1')
        PRIOR2m=input('According to Instructions 3, Section 3.4, retype the degrees of freedom for the covariance matrices at Time 2   ')
        [Ro,Co]=size(PRIOR2m);
    end
    while Co~=J(2)
        warning('The number of columns has to equal the number of clusters at time 2')
        PRIOR2m=input('According to Instructions 3, Section 3.4, retype the degrees of freedom for the covariance matrices at Time 2   ')
        [Ro,Co]=size(PRIOR2m);
    end
    
    %Time 3
    if Times==3
    PRIOR3m=input('According to Instruction 3, Section 3.4, specify the degrees of freedom for the covariance matrices at Time 3    ')
    %WARNINGS
    [Ro,Co]=size(PRIOR3m);
    while Ro~=1
        warning('This is a vector, i.e. the number of rows has to be equal to 1')
        PRIOR3m=input('According to Instructions 3, Section 3.4, retype the degrees of freedom for the covariance matrices at Time 3   ')
        [Ro,Co]=size(PRIOR3m);
    end
    while Co~=J(3)
        warning('The number of columns has to equal the number of clusters at time 3')
        PRIOR3m=input('According to Instructions 3, Section 3.4, retype the degrees of freedom for the covariance matrices at Time 3   ')
        [Ro,Co]=size(PRIOR3m);
    end
    end

    %Time 1
    for i=1:J(1)
        Covariance1=input(['According to Instructions 4, Section 3.4, specify prior covariance matrix at Time 1', num2str(i),' '])
        %WARNINGS    
        [Ro,Co]=size(Covariance1);
        while Ro~=D(1)
            warning('The number of rows has to equal the number of dimensions D(1)')
            Covariance1=input(['According to Instrcutions 4, Section 3.4, retype prior covariance matrix at Time 1 ', num2str(i),' '])
            [Ro,Co]=size(Covariance1);
        end
        while Co~=D(1)
            warning('The number of columns has to equal the number of dimensions D(1)')
            Covariance1=input(['According to Instrcutions 4, Section 3.4, retype prior covariance matrix at Time 1 ', num2str(i),' '])
            [Ro,Co]=size(Covariance1);
        end
        PRIOR1.sig(:,i*D(1)-(D(1)-1):i*D(1))=Covariance1;
        PRIOR1.psi(:,i*D(1)-(D(1)-1):i*D(1))=PRIOR1m(i)*Covariance1;
    end
    
    %Time 2
    for i=1:J(2)
        Covariance2=input(['According to Instructions 4, Section 3.4, specify prior covariance matrix at Time 2', num2str(i),' '])
        %WARNINGS    
        [Ro,Co]=size(Covariance2);
        while Ro~=D(2)
            warning('The number of rows has to equal the number of dimensions D(2)')
            Covariance2=input(['According to Instrcutions 4, Section 3.4, retype prior covariance matrix at Time 2 ', num2str(i),' '])
            [Ro,Co]=size(Covariance2);
        end
        while Co~=D(2)
            warning('The number of columns has to equal the number of dimensions D(2)')
            Covariance2=input(['According to Instrcutions 4, Section 3.4, retype prior covariance matrix at Time 2 ', num2str(i),' '])
            [Ro,Co]=size(Covariance2);
        end
        PRIOR2.sig(:,i*D(2)-(D(2)-1):i*D(2))=Covariance2;
        PRIOR2.psi(:,i*D(2)-(D(2)-1):i*D(2))=PRIOR2m(i)*Covariance2;
    end
    
    %Time 3
    if Times==3
    for i=1:J(3)
        Covariance3=input(['According to Instructions 4, Section 3.4, specify prior covariance matrix at Time 3', num2str(i),' '])
        %WARNINGS    
        [Ro,Co]=size(Covariance3);
        while Ro~=D(3)
            warning('The number of rows has to equal the number of dimensions D(3)')
            Covariance3=input(['According to Instrcutions 4, Section 3.4, retype prior covariance matrix at Time 3 ', num2str(i),' '])
            [Ro,Co]=size(Covariance3);
        end
        while Co~=D(3)
            warning('The number of columns has to equal the number of dimensions D(3)')
            Covariance3=input(['According to Instrcutions 4, Section 3.4, retype prior covariance matrix at Time 3 ', num2str(i),' '])
            [Ro,Co]=size(Covariance3);
        end
        PRIOR3.sig(:,i*D(3)-(D(3)-1):i*D(3))=Covariance3;
        PRIOR3.psi(:,i*D(3)-(D(3)-1):i*D(3))=PRIOR3m(i)*Covariance3;
    end
    end

    %PRIOR SPECIFICATIONS - CLUSTER PROBABILITIES AT TIME 1
    %Time 1
    PRIOR1alfa=input('According to Instructions 5, Section 3.4, specify the alfa parameters for the cluster proportions at Time 1')
    %WARNINGS    
    [Ro,Co]=size(PRIOR1alfa);
    while Ro~=1
        warning('This is a vector, i.e. the number of rows has to be equal to 1')
        PRIOR1alfa=input('According to Instructions 5, Section 3.4, retype the alfa parameters for the cluster proportions at Time 1 ')
        [Ro,Co]=size(PRIOR1alfa);
    end
    while Co~=J(1)
        warning('The number of columns has to equal the number of clusters at time 1')
        PRIOR1alfa=input('According to Instructions 5, Section 3.4, retype the alfa parameters for the cluster proportions at Time 1 ')
        [Ro,Co]=size(PRIOR1alfa);
    end
     
    %PRIOR SPECIFICATIONS - TRANSITION PROBABILITIES
    PRIOR1beta=input('According to Instructions 6, Section 3.4, specify the beta parameters for the transitions between Time 1 and 2  ')
       [Ro,Co]=size(PRIOR1beta);  
    while Ro~=J(1)
        warning('The number of rows has to equal the number of clusters at time 1')
        PRIOR1beta=input('According to Instructions 6, Section 3.4, retype the beta parameters for the transitions between Time 1 and 2 ')
        [Ro,Co]=size(PRIOR1beta);
    end
    while Co~=J(2)
        warning('The number of columns has to equal the number of clusters at time 2')
        PRIOR1beta=input('According to Instructions 6, Section 3.4, retype the beta parameters for the transitions between Time 1 and 2 ')
        [Ro,Co]=size(PRIOR1beta);
    end   
    
    if Times==3
        PRIOR2beta=input('According to Instructions 6, Section 3.4, specify the beta parameters for the transitions between Time 2 and 3  ')
       [Ro,Co]=size(PRIOR2beta);  
        while Ro~=J(2)
            warning('The number of rows has to equal the number of clusters at time 2')
            PRIOR2beta=input('According to Instructions 6, Section 3.4, retype the beta parameters for the transitions between Time 2 and 3 ')
            [Ro,Co]=size(PRIOR2beta);
        end
        while Co~=J(3)
            warning('The number of columns has to equal the number of clusters at time 3')
            PRIOR2beta=input('According to Instructions 6, Section 3.4, retype the beta parameters for the transitions between Time 2 and 3 ')
            [Ro,Co]=size(PRIOR2beta);
        end  
    end
end

PRIOR1.tau=PRIOR1tau;
PRIOR1.ksi=PRIOR1ksi;
PRIOR1.m=PRIOR1m;
PRIOR1.alfa=PRIOR1alfa;
PRIOR1.beta=PRIOR1beta;

PRIOR2.tau=PRIOR1tau;
PRIOR2.ksi=PRIOR1ksi;
PRIOR2.m=PRIOR1m;

if Times==3
PRIOR3.tau=PRIOR3tau;
PRIOR3.ksi=PRIOR3ksi;
PRIOR3.m=PRIOR3m;
PRIOR2.beta=PRIOR2beta;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%CREATES MATRICES WHICH WILL BE FILLED UP IN THE ITERATIONS
M1=zeros(D(1),J(1),T+1);
Theta1=zeros(T+1,J(1));
ystreck1=[];
f1=zeros(J(1),n);
Sigma1=zeros(D(1),J(1)*D(1),T+1);
Antal1=zeros(n,J(1));
V1=zeros(n);

M2=zeros(D(2),J(2),T+1);
Theta2=zeros(T+1,J(2));
ystreck2=[];
f2=zeros(J(2),n);
Sigma2=zeros(D(2),J(2)*D(2),T+1);
Antal2=zeros(n,J(2));
V2=zeros(n);
Vcount1=zeros(J(1),n);
Vcount2=zeros(J(2),n);
Q1=zeros(J(1),J(2),T+1);

%START VALUES GENERATED BY K-MEANS CLUSTERING
[IDX1, C1, SUMD1, Do1] = KMEANS(Y1',J(1));
V01=IDX1';

[IDX2, C2, SUMD2, Do2] = KMEANS(Y2',J(2));
V02=IDX2';

for j=1:J(1)           
    for d=1:D(1)       
        Theta01(1,j)=(length(find(V01==j)))/length(Y1);
        M01(d,j)=mean(Y1(d,V01==j));
        Sigma01(:,j*D(1)-(D(1)-1):j*D(1))=cov((Y1(:,V01==j))');
    end
end

for j=1:J(2)           
    for d=1:D(2)       
        Theta02(1,j)=(length(find(V02==j)))/length(Y2);
        M02(d,j)=mean(Y2(d,V02==j));
        Sigma02(:,j*D(2)-(D(2)-1):j*D(2))=cov((Y2(:,V02==j))');
    end
end

M1(:,:,1)=M01(:,1:J(1)); 
V1=V01;
Theta1(1,:)=Theta01;
Sigma1(:,:,1)=Sigma01(:,1:D(1)*J(1));

M2(:,:,1)=M02(:,1:J(2)); 
V2=V02;
Theta2(1,:)=Theta02;
Sigma2(:,:,1)=Sigma02(:,1:D(2)*J(2));

%Start values for transition matrices
Q1(:,:,1)=(1/J(2))*ones(J(1),J(2));

%START VALUES IF T=3
if Times==3
    M3=zeros(D(3),J(3),T+1);
    Theta3=zeros(T+1,J(3));
    ystreck3=[];
    f3=zeros(J(3),n);
    Sigma3=zeros(D(3),J(3)*D(3),T+1);
    Antal3=zeros(n,J(3));
    V3=zeros(n);
    Vcount3=zeros(J(3),n);
    Q2=zeros(J(2),J(3),T+1);

    [IDX3, C3, SUMD3, Do3] = KMEANS(Y3',J(3));
    V03=IDX3';
    
    for j=1:J(3)           
        for d=1:D(3)       
            Theta03(1,j)=(length(find(V03==j)))/length(Y3);
            M03(d,j)=mean(Y3(d,V03==j));
            Sigma03(:,j*D(3)-(D(3)-1):j*D(3))=cov((Y3(:,V03==j))');
        end
    end
    
    Q2(:,:,1)=(1/J(3))*ones(J(2),J(3));
    
    M3(:,:,1)=M03(:,1:J(3)); 
    V3=V03;
	Theta3(1,:)=Theta03;
    Sigma3(:,:,1)=Sigma03(:,1:D(3)*J(3));      
end

%THE SIMULATIONS START HERE
tic
for t=1:T
     
    for j=1:J(1)
        m1(j)=length(find(V1==j));
    end

    for j=1:J(2)
        m2(j)=length(find(V2==j));
    end

    if Times==3
    for j=1:J(3)
        m3(j)=length(find(V3==j));
    end
    end
    
    %SIMULATIONS OF COVARIANCE MATRICES
    %Cov 1
    df1=m1+PRIOR1.m;      % df för sigma är 10 och samma för alla j
    for j=1:J
    
        if m1(j)>0
            P1=Y1(:,V1==j);
            ystreck21=mean(P1,2);
            for i=1:m1(j)
                H1(:,:,i)=(P1(:,i)-ystreck21)*(P1(:,i)-ystreck21)';
            end
            W1=sum(H1,3);
            H1=[];
            P1=[];
            SigmaW1=PRIOR1.psi(:,j*D(1)-(D(1)-1):j*D(1))+W1+((m1(j)*PRIOR1.tau(j))/(m1(j)+PRIOR1.tau(j)))*(ystreck21-PRIOR1.ksi(:,j))*(ystreck21-PRIOR1.ksi(:,j))';
            temp1=(triu(SigmaW1,1))';
            SigmaW1=triu(SigmaW1)+temp1;
            Sigma1(:,j*D(1)-(D(1)-1):j*D(1),t+1)=iwishrnd((SigmaW1),df1(1,j));
        else  
            Sigma1(:,j*D(1)-(D(1)-1):j*D(1),t+1)=iwishrnd((PRIOR1.sig),PRIOR1.m(j));   
        end
    
    end

    %Cov 2
    df2=m2+PRIOR2.m;      % df för sigma är 10 och samma för alla j
    for j=1:J(2)
    
        if m2(j)>0
            P2=Y2(:,V2==j);
        	ystreck22=mean(P2,2);
            for i=1:m2(j)
                H2(:,:,i)=(P2(:,i)-ystreck22)*(P2(:,i)-ystreck22)';
            end 
            W2=sum(H2,3);
            H2=[];
            P2=[];
            SigmaW2=PRIOR2.psi(:,j*D(2)-(D(2)-1):j*D(2))+W2+((m2(j)*PRIOR2.tau(j))/(m2(j)+PRIOR2.tau(j)))*(ystreck22-PRIOR2.ksi(:,j))*(ystreck22-PRIOR2.ksi(:,j))';
            temp2=(triu(SigmaW2,1))';
            SigmaW2=triu(SigmaW2)+temp2;
            Sigma2(:,j*D(2)-(D(2)-1):j*D(2),t+1)=iwishrnd((SigmaW2),df2(1,j));
        else
            Sigma2(:,j*D(2)-(D(2)-1):j*D(2),t+1)=iwishrnd((PRIOR2.sig),PRIOR2.m(j));   
        end
    end

    %Cov 3
    if Times==3
    df3=m3+PRIOR3.m;      % df för sigma är 10 och samma för alla j
        for j=1:J(3)  
            if m3(j)>0
            P3=Y3(:,V3==j);
            ystreck23=mean(P3,2);
                for i=1:m3(j)
                    H3(:,:,i)=(P3(:,i)-ystreck23)*(P3(:,i)-ystreck23)';
                end
                W3=sum(H3,3);
                H3=[];
                P3=[];
                SigmaW3=PRIOR3.psi(:,j*D(3)-(D(3)-1):j*D(3))+W3+((m3(j)*PRIOR3.tau(j))/(m3(j)+PRIOR3.tau(j)))*(ystreck23-PRIOR3.ksi(:,j))*(ystreck23-PRIOR3.ksi(:,j))';
                temp3=(triu(SigmaW3,1))';
                SigmaW3=triu(SigmaW3)+temp3;
                Sigma3(:,j*D(3)-(D(3)-1):j*D(3),t+1)=iwishrnd((SigmaW3),df3(1,j));
            else
                Sigma3(:,j*D(3)-(D(3)-1):j*D(3),t+1)=iwishrnd((PRIOR3.sig),PRIOR3.m(j));   
            end
        end
    end

    %SIMULATION OF MEAN VECTORS
    
    %M1
    for j=1:J(1)
        if m1(j)>0
            ystreck1(:,j)=mean(Y1(:,V1==j),2);
            Mpost1(:,j)=(PRIOR1.tau(j)*PRIOR1.ksi(:,j)+(m1(j)*ystreck1(:,j)))/(PRIOR1.tau(j)+m1(j));
            Sigmapost1=Sigma1(:,j*D(1)-(D(1)-1):j*D(1),t+1)/(PRIOR1.tau(j)+m1(j));
            M1(:,j,t+1)=mvnrnd(Mpost1(:,j),Sigmapost1); 
        else
            M1(:,j,t+1)=mvnrnd(PRIOR1.ksi(:,j),PRIOR1.sig);
        end
    end

    %M2
    for j=1:J(2)
        if m2(j)>0
            ystreck2(:,j)=mean(Y2(:,V2==j),2);
            Mpost2(:,j)=(PRIOR2.tau(j)*PRIOR2.ksi(:,j)+(m2(j)*ystreck2(:,j)))/(PRIOR2.tau(j)+m2(j));
            Sigmapost2=Sigma2(:,j*D(2)-(D(2)-1):j*D(2),t+1)/(PRIOR2.tau(j)+m2(j));
            M2(:,j,t+1)=mvnrnd(Mpost2(:,j),Sigmapost2);
        else
            M2(:,j,t+1)=mvnrnd(PRIOR2.ksi(:,j),PRIOR2.sig);
        end
    end

    %M3
    if Times==3
        for j=1:J(3)
            if m3(j)>0
                ystreck3(:,j)=mean(Y3(:,V3==j),2);
                Mpost3(:,j)=(PRIOR3.tau(j)*PRIOR3.ksi(:,j)+(m3(j)*ystreck3(:,j)))/(PRIOR3.tau(j)+m3(j));
                Sigmapost3=Sigma3(:,j*D(3)-(D(3)-1):j*D(3),t+1)/(PRIOR3.tau(j)+m3(j));
                M3(:,j,t+1)=mvnrnd(Mpost3(:,j),Sigmapost3); 
            else
                M3(:,j,t+1)=mvnrnd(PRIOR3.ksi(:,j),PRIOR3.sig);
            end
        end
    end

    %SIMULATION OF CLUSTER PROBABILITIES         
    Postalfa1=PRIOR1.alfa+m1;
    G1=gamrnd(Postalfa1,1);        
    Theta1(t+1,:)=G1./sum(G1);
    %Theta1(t+1,:)=sort(Theta1(t+1,:),'descend');
    

    %SIMULATION OF THE CLASSIFICATION VECTORS V 
    for j=1:J(1)
        Inv1=inv(Sigma1(:,j*D(1)-(D(1)-1):j*D(1),t+1));
        for i=1:n 
            f1(j,i)=((2*pi)^(-D(1)/2)*(det(Sigma1(:,j*D(1)-(D(1)-1):j*D(1),t+1)))^(-1/2))*(exp(-1/2*(Y1(:,i)-M1(:,j,t+1))'*Inv1*(Y1(:,i)-M1(:,j,t+1))));
        end
    end      
    
    for j=1:J(2)
        Inv2=inv(Sigma2(:,j*D(2)-(D(2)-1):j*D(2),t+1));
        for i=1:n 
            f2(j,i)=((2*pi)^(-D(2)/2)*(det(Sigma2(:,j*D(2)-(D(2)-1):j*D(2),t+1)))^(-1/2))*(exp(-1/2*(Y2(:,i)-M2(:,j,t+1))'*Inv2*(Y2(:,i)-M2(:,j,t+1))));
        end
    end     
    
    if Times==3
        for j=1:J(3)
            Inv3=inv(Sigma3(:,j*D(3)-(D(3)-1):j*D(3),t+1));
            for i=1:n 
                f3(j,i)=((2*pi)^(-D(3)/2)*(det(Sigma3(:,j*D(3)-(D(3)-1):j*D(3),t+1)))^(-1/2))*(exp(-1/2*(Y3(:,i)-M3(:,j,t+1))'*Inv3*(Y3(:,i)-M3(:,j,t+1))));
            end
        end     
    end  
    
    if Times==2
        Number=zeros(J(1),J(2));
        for i=1:n
            for j=1:J(1)
                for k=1:J(2)
                    ksiTemp(j,k)=Theta1(t+1,j)*Q1(j,k,t)*f1(j,i)*f2(k,i);
                end
            end
     
        Summa=sum(sum(ksiTemp));

        ksi=ksiTemp./Summa;
        ksilong=reshape(ksi,1,J(1)*J(2));
        U=rand;
        Boris=U<cumsum(ksilong,2);
        VTemp=min(find(Boris));
        Boris=[];
        Jess=zeros(J(1),J(2));
        Jess(VTemp)=1;
        [Row,Col]=find(Jess==1);
        V1(i)=Row;
        V2(i)=Col;
    
        Number=Number+Jess;
        end
        
        if t>F
            for i=1:n
                Vcount1(V1(i),i)=Vcount1(V1(i),i)+1;
                Vcount2(V2(i),i)=Vcount2(V2(i),i)+1;
            end
        end
        
    else
  
        Number=zeros(J(1),J(2),J(3));
        for i=1:n
            if f1(1,1)*f2(1,1)*f3(1,1)==0
        
                for j=1:J(1)
                    for k=1:J(2)
                        for l=1:J(3)
                            ksiTemp(j,k,l)=Theta1(t+1,j)*Q1(j,k,t)*Q2(k,l,t)*f1(j,i)*10^100*f2(k,i)*10^100*f3(l,i)*10^100;  
                        end
                    end
                end
        
            else
        
                for j=1:J(1)
                    for k=1:J(2)
                        for l=1:J(3)
                            ksiTemp(j,k,l)=Theta1(t+1,j)*Q1(j,k,t)*Q2(k,l,t)*f1(j,i)*f2(k,i)*f3(l,i);
                        end
                    end
                end      
            end
         
            Test=all(ksiTemp);
            if Test==1
             
                for j=1:J(1)
                    for k=1:J(2)
                        for l=1:J(3)
                            ksiTemp(j,k,l)=10^10*Theta1(t+1,j)*Q1(j,k,t)*Q2(k,l,t)*f1(j,i)*f2(k,i)*f3(l,i);  
                        end
                    end
                end
            end
                              
            Summa=sum(sum(sum(ksiTemp)));

            ksi=ksiTemp./(Summa);
            ksilong=reshape(ksi,1,J(1)*J(2)*J(3));
            U=rand;
            Boris=U<cumsum(ksilong,2);
            VTemp=min(find(Boris));
            Boris=[];
            Jess=zeros(J(1),J(2),J(3));
            Jess(VTemp)=1; 
            B1=find(Jess==1);
        
            if B1<=J(1)*J(2)
                [Row,Col]=find(Jess==1);
                Lay=1;
              
            elseif J(1)*J(2)<B1 & B1<=(J(1)*J(2)*2)
                [Row,Col]=find(Jess(:,:,2)==1);
                Lay=2;
              
            else
                [Row,Col]=find(Jess(:,:,3)==1);
                Lay=3;
            end
        
            V1(i)=Row;
            V2(i)=Col;
            V3(i)=Lay;  
            Number=Number+Jess;
        end 
  
        if t>F
            for i=1:n
                Vcount1(V1(i),i)=Vcount1(V1(i),i)+1;
                Vcount2(V2(i),i)=Vcount2(V2(i),i)+1;
                Vcount3(V3(i),i)=Vcount3(V3(i),i)+1;
            end
        end
        
    end

    %SIMULATION OF TRANSITION MATRIX Q
    if Times==2
        for j=1:J(1)
            Postbeta1=PRIOR1.beta(j,:)+Number(j,:);
            Gam1=gamrnd(Postbeta1,1);        %Genererar J gammaförd. variabler med parametrar Postalfa
            Q1(j,:,t+1)=Gam1./sum(Gam1);
        end       
    else      
        Number1=sum(Number,3);
        for j=1:J(1)
            Postbeta1=PRIOR1.beta(j,:)+Number1(j,:);
            Gam1=gamrnd(Postbeta1,1);       
            Q1(j,:,t+1)=Gam1./sum(Gam1);
        end 
        Number2=sum(Number);
    
        for i=1:J(3)
            Number3(:,i)=Number2(:,:,i);
        end
    
        for j=1:J(2)
            Postbeta2=PRIOR2.beta(j,:)+Number3(j,:);
            Gam2=gamrnd(Postbeta2,1);       
            Q2(j,:,t+1)=Gam2./sum(Gam2);
        end
    end
    
    %CALCULATIONS OF CLUSTER PROBABILITIES
    Theta2(t+1,:)=Theta1(t+1,:)*Q1(:,:,t+1);
    if Times==3
        Theta3(t+1,:)=Theta2(t+1,:)*Q2(:,:,t+1); 
    end
    
end
toc

%SUMMARY STATISTICS
MEAN1=mean(M1(:,:,F+1:T),3)
MEAN2=mean(M2(:,:,F+1:T),3)
if Times==3
    MEAN3=mean(M3(:,:,F+1:T),3)
end

PROB1=mean(Theta1(F+1:T,:))
PROB2=mean(Theta2(F+1:T,:))
if Times==3
    PROB3=mean(Theta3(F+1:T,:))
end
          
Sigmakapad1=mean(Sigma1(:,:,F+1:T),3);
for j=1:J(1)
    eval(['COV1' num2str(j) ' = Sigmakapad1(:,j*D(1)-(D(1)-1):j*D(1))']);
end

Sigmakapad2=mean(Sigma2(:,:,F+1:T),3);
for j=1:J(2)
    eval(['COV2' num2str(j) ' = Sigmakapad2(:,j*D(2)-(D(2)-1):j*D(2))']);
end

if Times==3
    Sigmakapad3=mean(Sigma3(:,:,F+1:T),3);
    for j=1:J(3)
        eval(['COV3' num2str(j) ' = Sigmakapad3(:,j*D(3)-(D(3)-1):j*D(3))']);
    end
end

TRANSQ1=mean(Q1(:,:,F+1:T),3)
if Times==3
    TRANSQ2=mean(Q2(:,:,F+1:T),3)
end

for i=1:n
    for j=1:J(1)
        CLUSTERPROB1(i,j)=(Vcount1(j,i))/(T-F);
    end
end  

for i=1:n
    for j=1:J(2)
        CLUSTERPROB2(i,j)=(Vcount2(j,i))/(T-F);
    end
end  

if Times==3
    for i=1:n
        for j=1:J(3)
            CLUSTERPROB(i,j)=(Vcount3(j,i))/(T-F);
        end
    end
end

