%MODEL SPECIFICATIONS  
[K,n] = size(Y);     %K is the dimension of data and n is the number of observations 
if K>n
    warning('The number of observations n have to exceed the dimension K of data. Respecify the data matrix Y')
    break
end
    
J=input('Give the number of clusters J ')
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  ')

if Default==0     
    PRIORtau=ones(1,J);
    PRIORksi=repmat(mean(Y')',1,J);
    PRIORm=10*ones(1,J);
    Covariance=cov(Y');
    for i=1:J
        PRIOR.sig(:,i*K-(K-1):i*K)=Covariance;
        PRIOR.psi(:,i*K-(K-1):i*K)=PRIORm(i)*Covariance;
    end
    PRIORalfa=5*ones(1,J);
else

    %PRIOR SPECIFICATIONS - CLUSTER MEANS
    %Specify the precision parameters in vector form. The length of the vector is equal to J
    PRIORtau=input('According to Instuction 1, Section 3.1, specify the precision parameters for the mean variables   ') 
    %WARNINGS
    [Ro,Co]=size(PRIORtau);
    while Ro~=1
        warning('This is a vector, i.e. the number of rows has to be equal to 1')
        PRIORtau=input('According to Instuction 1, Section 3.1, retype the precision parameters for the mean variables  ') 
        [Ro,Co]=size(PRIORtau);
    end
    while Co~=J
        warning('The number of columns has to equal the number of clusters J')
        PRIORtau=input('According to Instuction 1, Section 3.1, retype the precision parameters for the mean variables  ') 
        [Ro,Co]=size(PRIORtau);
    end

    %Give the prior specification of the mean values for each cluster in all dimensions in matrix form, 
    %Rows=dimensions=K and columns=clusters=J
    PRIORksi=input('According to Intruction 2, Section 3.2, specify the prior of the mean values for all clusters   ')
    %WARNINGS
    [Ro,Co]=size(PRIORksi);
    while Ro~=K
        warning('The number of rows has to equal the number of variables K')
        PRIORksi=input('According to Intruction 2, Section 3.2, retype the prior of the mean values for all clusters  ')
        [Ro,Co]=size(PRIORksi);
    end
    while Co~=J
        warning('The number of columns has to equal the number of clusters J')
        PRIORksi=input('According to Intruction 2, Section 3.2, retype the prior of the mean values for all clusters  ')
        [Ro,Co]=size(PRIORksi);
    end

    %PRIOR SPECIFICATIONS - CLUSTER COVARIANCES
    %Give the prior specification of mj, the degrees of freedom in vector form. The length of the vector is equal to J.
    PRIORm=input('According to Instructions 3, Section 3.3, specify the degrees of freedom for the covariance matrices    ')
    %WARNINGS
    [Ro,Co]=size(PRIORm);
    while Ro~=1
        warning('This is a vector, i.e. the number of rows has to be equal to 1')
        PRIORm=input('According to Instructions 3, Section 3.3, retype the degrees of freedom for the covariance matrices   ')
        [Ro,Co]=size(PRIORm);
    end
    while Co~=J
        warning('The number of columns has to equal the number of clusters J')
        PRIORm=input('According to Instructions 3, Section 3.3, retype the degrees of freedom for the covariance matrices   ')
        [Ro,Co]=size(PRIORm);
    end

    %Give the prior specification of the sigma matrices. One matrix for each cluster J
    for i=1:J
        Covariance=input(['According to Instructions 4, Section 3.4, specify prior covariance matrix ', num2str(i),' '])
        %WARNINGS    
        [Ro,Co]=size(Covariance);
        while Ro~=K
            warning('The number of rows has to equal the number of dimensions K')
            Covariance=input(['According to Instrcutions 4, Section 3.4, retype prior covariance matrix ', num2str(i),' '])
            [Ro,Co]=size(Covariance);
        end
        while Co~=K
            warning('The number of columns has to equal the number of dimensions K')
            Covariance=input(['According to Instrcutions 4, Section 3.4, retype prior covariance matrix ', num2str(i),' '])
            [Ro,Co]=size(Covariance);
        end
        PRIOR.sig(:,i*K-(K-1):i*K)=Covariance;
        PRIOR.psi(:,i*K-(K-1):i*K)=PRIORm(i)*Covariance;
    end

    %PRIOR SPECIFICATIONS - CLUSTER PROBABILITIES
    %Give the prior specification of the alfa parameters in vector form. The length of the vector is equal to J.
    PRIORalfa=input('According to Instructions 5, Section 3.5, specify the alfa parameters for the cluster proportions ')
    %WARNINGS    
    [Ro,Co]=size(PRIORalfa);
    while Ro~=1
        warning('This is a vector, i.e. the number of rows has to be equal to 1')
        PRIORalfa=input('According to Instructions 5, Section 3.5, retype the alfa parameters for the cluster proportions ')
        [Ro,Co]=size(PRIORalfa);
    end
    while Co~=J
        warning('The number of columns has to equal the number of clusters J')
        PRIORalfa=input('According to Instructions 5, Section 3.5, retype the alfa parameters for the cluster proportions ')
        [Ro,Co]=size(PRIORalfa);
    end
end

PRIOR.tau=PRIORtau;
PRIOR.ksi=PRIORksi;
PRIOR.m=PRIORm;
PRIOR.alfa=PRIORalfa;

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

%CREATES MATRICES WHICH WILL BE FILLED UP IN THE ITERATIONS
M=zeros(K,J,T+1);
Theta=zeros(T+1,J);
f=zeros(J,n);
Sigma=zeros(K,J*K,T+1);
Vcount=zeros(J,n);

%START VALUES GENERATED BY K-MEANS CLUSTERING
[IDX, C, SUMD, D] = KMEANS(Y', J);
V0=IDX';

for j=1:J           
    for k=1:K       
        Theta0(1,j)=(length(find(V0==j)))/length(Y);
        M0(k,j)=mean(Y(k,V0==j));
        Sigma0(:,j*K-(K-1):j*K)=cov((Y(:,V0==j))');
    end
end

M(:,:,1)=M0(:,1:J); 
V(1,:)=V0;
Theta(1,:)=Theta0;
Sigma(:,:,1)=Sigma0(:,1:K*J);

%THE SIMULATIONS START HERE
tic
for t=1:T   
    
    for j=1:J
        m(j)=length(find(V==j));   %The number of observations in each cluster
    end
 
    %SIMULATIONS OF COVARIANCE MATRICES
    df=m+PRIOR.m;     
    for j=1:J
        if  m(j)>0                      %to be able to count the cluster means there need to be at leat 1 observation in each cluster
            P=Y(:,V==j);
            ystreck2=mean(P,2);

            for i=1:m(j)
                H(:,:,i)=(P(:,i)-ystreck2)*(P(:,i)-ystreck2)';
            end
            
            W=sum(H,3);
            H=[];
            P=[];

            SigmaW=PRIOR.psi(:,j*K-(K-1):j*K)+W+((m(j)*PRIOR.tau(j))/(m(j)+PRIOR.tau(j)))*(ystreck2-PRIOR.ksi(:,j))*(ystreck2-PRIOR.ksi(:,j))';
            temp=(triu(SigmaW,1))';
            SigmaW=triu(SigmaW)+temp;
            Sigma(:,j*K-(K-1):j*K,t+1)=iwishrnd((SigmaW),df(1,j));
        else                            %if one cluster is empty the conditional distribution is the same as the prior
            Sigma(:,j*K-(K-1):j*K,t+1)=iwishrnd(PRIOR.sig(:,j*K-(K-1):j*K),PRIOR.m(j));  
        end
    end

    %SIMULATION OF MEAN VECTORS
    for j=1:J
        if  m(j)>0                      %to be able to count the cluster means there need to be at least 1 observation in each cluster
            ystreck(:,j)=mean(Y(:,V==j),2);
            Mpost(:,j)=(PRIOR.tau(j)*PRIOR.ksi(:,j)+(m(j)*ystreck(:,j)))/(PRIOR.tau(j)+m(j));
            Sigmapost=Sigma(:,j*K-(K-1):j*K,t+1)/(PRIOR.tau(j)+m(j));
            M(:,j,t+1)=mvnrnd(Mpost(:,j),Sigmapost);  
        else
            M(:,j,t+1)=mvnrnd(PRIOR.ksi(:,j),PRIOR.sig(:,j*K-(K-1):j*K));  %if one cluster is empty the conditional distribution is the same as the prior
        end
    end

    %SIMULATION OF CLUSTER PROBABILITIES            
    Postalfa=PRIOR.alfa+m;
    G=gamrnd(Postalfa,1);               
    Theta(t+1,:)=G./sum(G);
    %Theta(t+1,:)=sort(Theta(t+1,:),'descend');

    %SIMULATION OF THE CLASSIFICATION VECTOR V  
    for j=1:J
        Inv=inv(Sigma(:,j*K-(K-1):j*K,t+1));
        for i=1:n 
            f(j,i)=((2*pi)^(-K/2)*(det(Sigma(:,j*K-(K-1):j*K,t+1)))^(-1/2))*(exp(-1/2*(Y(:,i)-M(:,j,t+1))'*Inv*(Y(:,i)-M(:,j,t+1))));
        end
    end             
         
    Summa=Theta(t+1,:)*f;

    for i=1:n
        for j=1:J
            TT(i,j)=(Theta(t+1,j)*f(j,i))/Summa(i);
        end
    end

    U=rand(n,1);
    for i=1:n
        Boris=U(i)<cumsum((TT(i,:)),2);
        V(i)=min(find(Boris));
        Boris=[];
    end
    
    %Saves cluster classifications for all individuals after F iterations (burn-in)
    if t>F
        for i=1:n
            Vcount(V(i),i)=Vcount(V(i),i)+1;
        end
    end
   
end
toc

%SUMMARY STATISTICS
MEAN=mean(M(:,:,F+1:T),3)

PROB=mean(Theta(F+1:T,:))

SIGMA=mean(Sigma(:,:,F+1:T),3);
for j=1:J
eval(['Cov' num2str(j) ' = SIGMA(:,j*K-(K-1):j*K)'])
end

%Cluster probabilities for all observations
for i=1:n
    for j=1:J
        CLUSTERPROB(i,j)=(Vcount(j,i))/(T-F);
    end
end


