% Copyright (c) 2002, Daniel B. Rowe
% This Matlab program performs the computations 
% for a Multivariate Regression.
% It reads in the data file 'regdata.txt' and 
% computes estimates the parameters as in Rowe 2002.
% There are references listed at the end.

clear

% Have data arranged so that n is the number 
% of rows and p+q+1 is the number of columns
% Each row is an observation vector.

n=32;
p=3;
q=2;
 
disp(' ')
disp('Reading in the data matrix')
load 'regdata.txt'    % Reads in data from regdata.txt
X=regdata(:,1:p)
U=regdata(:,p+1:p+q+1)
clear regdata

% Computing estimated regression coefficients, MLE
disp(' ')
disp('Computing estimated regression coefficients')  
Bhat=X'*U*inv(U'*U)

% Computing estimated error covariance matrix, MLE 
disp(' ')
disp('Computing estimated error covariance matrix')
G=(X-U*Bhat')'*(X-U*Bhat')
Sigmahat=G/n


W=zeros(q+1,q+1);
W=inv(U'*U);
disp(' ')
disp('Computing marginal posterior statistics')
% Computing marginal t-statistics for each coefficient
t=zeros(p,q+1);
for col=1:q+1
   for row=1:p
    t(row,col)=Bhat(row,col)/sqrt(W(col,col)*G(row,row)/(n-q-1)) ; 
   end
end
t

%save reg


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% References                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Rowe, D.B. (2003). Multivariate Bayesian Statistics: 
% Models for Source Separation and Signal Unmixing.
% CRC Press, Boca Raton, Florida.

% Rowe, D.B. and Morgan, S.W. (2002). Computing FMRI 
% Activations: Coefficients and t-Statistics by Detrending 
% and Multiple Regression. Technical Report 39, Division 
% of Biostatistics, Medical College of Wisconsin, Milwaukee 
% Wisconsin, 53226. 

% Press, S. J. (1982). Applied Multivariate Analysis: Using 
% Bayesian and Frequentist Methods of Inference. Robert E. Krieger
% Publishing Company, Malabar, Florida. 

   12.1814    9.3538   20.9415    1    1    1
   15.8529   11.7987   28.0296    1    2    1
   20.5458   15.1322   35.0787    1    3    1
   23.9659   18.0548   42.3609    1    4    1
   28.0285   20.7695   48.9123    1    5    1
   32.2667   23.4573   56.1558    1    6    1
   36.0148   26.9852   63.1998    1    7    1
   39.9761   29.7473   70.2352    1    8    1
   31.7919   23.1536   58.7520    1    9   -1
   36.0736   26.1269   66.0530    1   10   -1
   39.6660   29.4231   73.0595    1   11   -1
   44.1786   32.1478   79.7481    1   12   -1
   48.4059   34.8391   86.8145    1   13   -1
   51.8271   38.0951   94.2706    1   14   -1
   56.2145   40.7477  100.9671    1   15   -1
   60.3135   43.9951  108.0975    1   16   -1
   75.6016   56.9879  133.0220    1   17    1
   79.6398   60.0000  139.8411    1   18    1
   84.1428   62.9205  146.8601    1   19    1
   87.9000   66.2738  154.1109    1   20    1
   92.1725   68.5315  160.7625    1   21    1
   96.2039   72.1070  168.1953    1   22    1
  100.1780   75.2239  175.1422    1   23    1
  104.3226   78.1827  181.7946    1   24    1
   96.1672   71.1445  170.9336    1   25   -1
  100.2977   74.0101  177.7031    1   26   -1
  103.6994   77.1693  184.4494    1   27   -1
  107.9951   80.1422  192.2466    1   28   -1
  111.9608   82.9361  198.8703    1   29   -1
  115.5990   85.9056  206.0818    1   30   -1
  120.0643   88.9260  213.0585    1   31   -1
  123.7359   91.6312  220.0054    1   32   -1

% Copyright (c) 2002, Daniel B. Rowe
% This Matlab program performs the computations 
% for a Bayesian Multivariate Regression.
% Conjugate prior distributions are used.
% This program will work on the student version.
% It reads in the data file 'sugar.txt' and 
% computes estimates the parameters as in Rowe 2002.
% Data is available to assess the hyperparameters from,
% the file which contains it is sugar0.txt.
% There are references listed at the end.

clear

n=36;
p=6;
q=3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hyperparameter assessment                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
assess='previus';
%assess='subject';

%%%%
if assess=='subject',
%%%%
disp('Subjectively Assessed Hyperparameters')

B0=[
  101.0000,    0.0556,   -0.3889,    0.8889;
   58.6944,    0.2500,    0.5833,    1.0278;
   66.3889,    2.5556,   -2.8889,    0.6111;
   76.9444,    3.0556,   -6.6111,   -0.9444;
   95.5278,    2.6944,   -6.3056,    1.5278;
  100.2222,    2.6111,   -2.5556,    2.7222]
D=[0.0278,          0,         0,         0;
        0,    0.0278,          0,         0;
        0,          0,    0.0278,         0;
        0,          0,         0,    0.0278]
Q=[ 109.20,     29.28,    -21.31,     10.26,    69.04,      66.89;
     29.28,    233.31,    201.53,    200.37,    107.54,     81.99;
    -21.31,    201.53,    289.36,    269.59,    107.33,     74.90;
     10.26,    200.37,    269.59,    461.22,    297.59,    205.36;
     69.04,    107.54,    107.33,    297.59,    478.62,    271.07;
     66.89,    081.99,     74.90,    205.36,    271.07,    248.49]
nu=36

%%%%
else
%%%%  
disp('Assessing hyperparameters from previous data.')
   
% Reads in data from sugar0.txt                 
% Have data arranged so that n0 is the number 
% of rows and p+q+1 is the number of columns
% and is of the form [X,U].
% Each row is an observation vector.

disp(' ')
disp('Reading in the previous data matrix')
load 'sugar0.txt';     % Reads in the data from sugar0.txt

disp(' ')
disp('The data matrix X0 is') 
X0=sugar0(:,1:p);
sizeX0=size(X0);
nu=sizeX0(1,1)

p=sizeX0(1,2)

disp(' ')
disp('The design matrix U is')
U0=sugar0(:,p+1:p+q+1);
sizeU0=size(U0);
nu=sizeU0(1,1)
clear sugar0

% Assessing prior hyperparameters B0 and D 
% for regression coefficient matrix 
disp(' ')
disp('Assessing prior mean B0 for regression coefficients')  
B0=X0'*U0*inv(U0'*U0)
D=inv(U0'*U0)

% Assessing prior hyperparameter matrix Q 
% for error covariance matrix 
disp(' ')
disp('Assessing prior scale matrix Q for error covariance matrix')
Q=(X0-U0*B0')'*(X0-U0*B0')

%%%%
end
%%%%

Dinv=zeros(q+1,q+1);
Dinv=inv(D);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computing marginal posterior mean estimates   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Computing marginal posterior mean estimates')

% Reads in data from bregdata.txt 
disp(' ')
disp('Reading in the data')
load 'sugar.txt'   % Reads in the data from sugar.txt

disp(' ')
disp('The data matrix X is') 
X=sugar(:,1:p)


disp(' ')
disp('The design matrix U is')
U=sugar(:,p+1:p+q+1)
clear sugar

disp(' ')
disp('Computing marginal posterior mean estimates')
% Computing marginal mean posterior estimates  
% for the regression coefficient matrix 

W=zeros(q+1,q+1);
W=inv(Dinv+U'*U);
Bbar=(X'*U+B0*Dinv)*W

% Computing marginal mean posterior estimates 
% for the error covariance matrix

G=Q+X'*X+B0*Dinv*B0'-(X'*U+B0*Dinv)*W*(X'*U+B0*Dinv)';
Sigmabar=G/(n+nu-2*p-2)

disp(' ')
disp('Computing marginal posterior statistics')
% Computing marginal t-statistics for each coefficient
t=zeros(p,q+1);
for col=1:q+1
   for row=1:p
    t(row,col)=Bbar(row,col)/sqrt(W(col,col)*G(row,row)/(n+nu-q-2*p+1)) ; 
   end
end
t


disp(' ')
disp('Computing maximum a posteriori estimates')
% Computing maximum a posteriori estimate  
% for the regression coefficient matrix 
  
Btilde=Bbar

% Computing maximum a posteriori estimate 
% for the error covariance matrix 
Sigmatilde=G/(n+nu+q+1)


disp(' ')
disp('Computing marginal conditional posterior statistics')
% Computing marginal z-statistics for each coefficient
z=zeros(p,q+1);
for col=1:q+1
   for row=1:p
    z(row,col)=Btilde(row,col)/sqrt(W(col,col)*G(row,row)/(n+nu+q+1)) ; 
   end
end
z


%save bregnat


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% References                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Rowe, D.B. (2003). Multivariate Bayesian Statistics: 
% Models for Source Separation and Signal Unmixing.
% CRC Press, Boca Raton, Florida.

% Rowe, D.B. (2001). A Bayesian Unobservable/Observable 
% Source Separation Model and Activation Determination in 
% FMRI. Social Science Working Paper 1120, Division of 
% Humanities and Social Sciences, California Institute 
% of Technology, Pasadena, CA 91125.

% Press S. J. (1989). Bayesian Statistics: Principles Models, 
% and Applications, John Wiley and Sons, New York.

% Press, S. J. (1982). Applied Multivariate Analysis: Using 
% Bayesian and Frequentist Methods of Inference. Robert E. Krieger
% Publishing Company, Malabar, Florida.

    96    54    61    63    93   103     1     1     1     1
    98    57    63    75    99   104     1     1     1     1
   104    77    88    91   113   110     1     1     1     1
   109    63    60    67    85   109     1     1     1     1
    98    59    65    72    95   103     1     1     1     1
   104    59    62    74    89    97     1     1     1     1
    97    63    70    72   101   102     1     1     1     1
   101    54    64    77    97   100     1     1     1     1
   107    59    67    61    69    99     1     1     1     1
    96    63    81    97   101    97     1     1    -1    -1
    99    48    70    94   108   104     1     1    -1    -1
   102    61    78    81    99   104     1     1    -1    -1
   112    67    76   100   112   112     1     1    -1    -1
    92    49    59    83   104   103     1     1    -1    -1
   101    53    63    86   104   102     1     1    -1    -1
   105    63    77    94   111   107     1     1    -1    -1
    99    61    74    76    89    92     1     1    -1    -1
    99    51    63    77    99   103     1     1    -1    -1
    98    53    62    71    81   101     1    -1     1    -1
   103    62    65    96   101   105     1    -1     1    -1
   102    54    60    57    64    69     1    -1     1    -1
   108    83    67    80   106   108     1    -1     1    -1
    92    56    60    61    73    79     1    -1     1    -1
   102    61    59    71    91   101     1    -1     1    -1
    94    51    53    55    86    83     1    -1     1    -1
    95    55    58    59    71    85     1    -1     1    -1
   103    47    59    64    92   100     1    -1     1    -1
   120    46    44    58   118   108     1    -1    -1     1
    95    65    75    85    96    95     1    -1    -1     1
    99    59    73    82   109   109     1    -1    -1     1
   105    50    58    84   107   107     1    -1    -1     1
    97    67    89   104   118   118     1    -1    -1     1
    97    46    50    59    78    91     1    -1    -1     1
   102    63    67    74    83    98     1    -1    -1     1
   104    69    81    98   104   105     1    -1    -1     1
   101    65    69    72    93    95     1    -1    -1     1

    96    37    31    33    35    41     1     1     1     1
    90    47    48    55    68    89     1     1     1     1
    99    49    55    64    74    97     1     1     1     1
    95    33    37    43    63    92     1     1     1     1
   107    62    62    85   110   117     1     1     1     1
    81    40    43    45    49    55     1     1     1     1
    95    49    56    63    68    88     1     1     1     1
   105    53    57    69   103   106     1     1     1     1
    97    50    53    59    82    96     1     1     1     1
    97    54    57    66    80    89     1     1    -1    -1
   105    66    83    95    97   100     1     1    -1    -1
   105    49    54    56    70    90     1     1    -1    -1
   106    79    92    95    99   100     1     1    -1    -1
    92    46    51    57    73    91     1     1    -1    -1
    91    61    64    71    80    90     1     1    -1    -1
   101    51    63    91    95    96     1     1    -1    -1
    87    53    55    57    78    89     1     1    -1    -1
    94    57    70    81    94    96     1     1    -1    -1
    98    48    55    71    91    96     1    -1     1    -1
    98    41    43    61    89   101     1    -1     1    -1
   103    60    56    61    76    97     1    -1     1    -1
    99    36    43    57    89   102     1    -1     1    -1
    97    44    51    58    85   105     1    -1     1    -1
    95    41    45    49    59    78     1    -1     1    -1
   109    65    62    72    93   104     1    -1     1    -1
    91    57    60    61    67    83     1    -1     1    -1
    99    43    48    52    61    86     1    -1     1    -1
   102    51    56    81    97   103     1    -1    -1     1
    96    57    55    72    85    89     1    -1    -1     1
   111    84    83    91   101   102     1    -1    -1     1
   105    57    67    83   100   103     1    -1    -1     1
   105    57    61    70    90    98     1    -1    -1     1
    98    55    67    88    94    95     1    -1    -1     1
    98    69    72    89    98    98     1    -1    -1     1
    90    53    61    78    94    95     1    -1    -1     1
   100    60    63    67    77   104     1    -1    -1     1

% Copyright (c) 2002, Daniel B. Rowe
% This Matlab program performs the computations 
% for Bayesian Factor Analysis as described
% in Rowe 2002 with conjugate priors.
% This program will work on the student version.
% It reads in the data file 'applicant.txt' and 
% with the hyperparameters in this program 
% estimates the parameters.
% There are references listed at the end.

clear

% Have data arranged so that n is the number 
% of rows and p is the number of columns
% Each row is an observation vector.

load 'applicant.txt'    % Reads in data from applicant.txt
[n,p]=size(applicant);
X=applicant;
clear applicant
xdsigma=sqrt(diag(cov(X))); % needed to scale data

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hyperparameter assessment                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Subjectively Assessed Hyperparameters')

mu0=zeros(p,1);
mu0=7.5*ones(p,1);
m=4;
Lambda0=zeros(p,m); % Lambda0 is a correlation matrix
Lambda0=.7 * [
0, 0, 1, 0;
0, 0, 0, 0;
0, 1, 0, 0;
0, 0, 0, 1;
1, 0, 0, 0;
1, 0, 0, 0;
0, 0, 0, 1;
1, 0, 0, 0;
0, 0, 1, 0;
1, 0, 0, 0;
1, 0, 0, 0;
1, 0, 0, 0;
1, 0, 0, 0;
0, 0, 0, 0;
0, 0, 1, 0];
C0=zeros(p,m);
C0=[mu0,Lambda0];
C0(:,1)=C0(:,1)./xdsigma; % Scale the prior mean
D0=1/10;
D=D0*eye(m+1);
nu=33;     % nu must be a positive integer
q0=0.200;
Q=q0*eye(p);


% These numbers can be changed depending on the data set.
ICMiterations=25000;      % # of ICM iterations
Gibbsburn=10000;          % # of burn in variates
Gibbsiterations=25000;    % # of samples for estimates


% Do not change any thing below here.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimation                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dinv=zeros(m+1,m+1);
Dinv=inv(D);

% Both the ICM and Gibbs sampling techniques need
% initial starting values. 
% Press/Shigemasu 1989 BFA estimates of parameters
muPS89=mean(X);
FPS89=zeros(n,m);
XmuPS89=X-ones(n,1)*muPS89;
XmuPS89=XmuPS89./(ones(n,1)*xdsigma');
W=XmuPS89'*XmuPS89+Q+Lambda0*inv(D(2:m+1,2:m+1))*Lambda0';
FPS891=(eye(n)-XmuPS89*inv(XmuPS89'*XmuPS89-W)*XmuPS89');
FPS89=FPS891*XmuPS89*inv(W)*Lambda0*inv(D(2:m+1,2:m+1))
LambdaPS89=(XmuPS89'*FPS89+Lambda0*inv(D(2:m+1,2:m+1)))*inv(inv(D(2:m+1,2:m+1))+FPS89'*FPS89)
PsiPS89mean=((XmuPS89-FPS89*LambdaPS89')'*(XmuPS89-FPS89*LambdaPS89')+(LambdaPS89-Lambda0)*inv(D(2:m+1,2:m+1))*(LambdaPS89-Lambda0)'+Q)/(n+m+nu-2*p-2)
PsiPS89mode=PsiPS89mean*(n+m+nu-2*p-2)/(n+m+nu)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ICM algorithm estimates                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize variables
FICM=zeros(n,m); 
ZICM=zeros(n,m+1);
CICM=zeros(p,m+1);
SigmaICM=zeros(p,p);

% Initial value
FICM=FPS89;

% Scale X matrix
X=X./(ones(n,1)*xdsigma');

for count=1:ICMiterations  
 ZICM=[ones(n,1),FICM];
 CICM= (X'*ZICM+C0*Dinv)*inv(Dinv+ZICM'*ZICM);
 SigmaICM=[(X-ZICM*CICM')'*(X-ZICM*CICM')+(CICM-C0)*Dinv*(CICM-C0)'+Q]/(n+nu+m+1);
 FICM1=(X-ones(n,1)*CICM(:,1)')*inv(SigmaICM)*CICM(:,2:m+1);
 FICM=FICM1*inv(eye(m)+CICM(:,2:m+1)'*inv(SigmaICM)*CICM(:,2:m+1)); 
end
 CICM(:,1)=CICM(:,1).*xdsigma; % Rescale the posterior mean
 CICM                          % Posterior estimate of Lambda0 is a matrix of correlations
 SigmaICM
 FICM
 
%%%ICM Significance
ZICM=[ones(n,1),FICM];
CICM(:,1)=CICM(:,1)./xdsigma; % Scale the posterior mean
WICM=zeros(m+1,m+1);
WICM=inv(Dinv+ZICM'*ZICM);

disp(' ')
disp('Computing marginal conditional posterior statistics')
% Computing marginal t-statistics for each coefficient
zICM=zeros(p,m+1);
for col=1:m+1
   for row=1:p
    zICM(row,col)=CICM(row,col)/sqrt(WICM(col,col)*SigmaICM(row,row)) ; 
   end
end
zICM
CICM(:,1)=CICM(:,1).*xdsigma; % Rescale the posterior mean


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gibbs sampling algorithm estimates             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize variables 
% The starting values don't matter
% But convergence is faster with good values
FGibbs=FICM;
CGibbs=CICM;
SigmaGibbs=SigmaICM;
ZGibbs=[ones(n,1),FGibbs];

Ftot=zeros(n,m);
Ctot=zeros(p,m+1);
C2tot=zeros(p*(m+1),p*(m+1));
Sigmatot=zeros(p,p);

YC = zeros(p,m+1);
YP = zeros(n+nu+m+1+p+1,p);
FGibbs=zeros(n,m);
WGibbs=zeros(m+1,m+1);
WGibbs=inv(Dinv+ZGibbs'*ZGibbs);

for count=1:Gibbsburn+Gibbsiterations   
  MC=(C0*Dinv+X'*ZGibbs)*WGibbs;
  YC = randn(p,m+1);
  AtC=chol(SigmaGibbs); 
  BtC=chol(WGibbs);
  CGibbs=AtC'*YC*BtC+MC;
   
  YP = randn(n+nu+m+1+p+1,p);
  AtP=chol((X-ZGibbs*CGibbs')'*(X-ZGibbs*CGibbs')+(CGibbs-C0)*Dinv*(CGibbs-C0)'+Q);
  SigmaGibbs=AtP'*inv(YP'*YP)*AtP;
  
  MF1=(X-ones(n,1)*CGibbs(:,1)')*inv(SigmaGibbs)*CGibbs(:,2:m+1);
  MF=MF1*inv(eye(m)+CGibbs(:,2:m+1)'*inv(SigmaGibbs)*CGibbs(:,2:m+1));
  FGibbs=randn(n,m);
  BtF=chol(inv(eye(m)+CGibbs(:,2:m+1)'*inv(SigmaGibbs)*CGibbs(:,2:m+1)));
  FGibbs=FGibbs*BtF+MF;
  
  ZGibbs=[ones(n,1),FGibbs];
  WGibbs=inv(Dinv+ZGibbs'*ZGibbs);
  
  if count>Gibbsburn
   Ftot=Ftot+FGibbs;
   Ctot=Ctot+CGibbs;
   C2tot=C2tot+vec(CGibbs)*vec(CGibbs)';
   Sigmatot=Sigmatot+SigmaGibbs;
  end
  
end


FGibbs=Ftot/Gibbsiterations
CGibbs=Ctot/Gibbsiterations;
DeltaGibbs=C2tot/Gibbsiterations-vec(CGibbs)*vec(CGibbs)';
CGibbs(:,1)=CGibbs(:,1).*xdsigma; % unscale the mean
CGibbs
SigmaGibbs=Sigmatot/Gibbsiterations 

%%%Gibbs Significance
disp(' ')
disp('Computing marginal posterior statistics')
% Computing marginal z-statistics for each coefficient
zGibbs=zeros(p,m+1);
for col=1:m+1
   for row=1:p
    zGibbs(row,col)=CGibbs(row,col)/sqrt(DeltaGibbs(p*(col-1)+row,p*(col-1)+row)); 
   end
end
zGibbs

%save bfanat


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% References                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Rowe, D.B. (2003). Multivariate Bayesian Statistics: 
% Models for Source Separation and Signal Unmixing.
% CRC Press, Boca Raton, Florida.

% Rowe, D. B. (2003). On Using the Sample Mean in Bayesian 
% Factor Analysis. Forthcoming in the Journal of 
% Interdisciplinary Mathematics, Vol. 6, No. 3, pp. 319-329. 

% Rowe, D.B. (2001). A Model for Bayesian Factor 
% Analysis with Jointly Distributed Means and Loadings. 
% Social Science Working Paper 1108, Division of Humanities 
% and Social Sciences, California Institute of Technology, 
% Pasadena, CA 91125. 

% Rowe, D. B. (2000). Incorporating Prior Knowlege Regarding 
% the Mean in Bayesian Factor Analysis. Social Science Working 
% Paper 1097, Division of Humanities and Social Sciences, 
% California Institute of Technology, Pasadena, CA 91125. 

% Rowe, D. B. (2000). On Estimating the Mean in Bayesian 
% Factor Analysis. Social Science Working Paper 1096, 
% Division of Humanities and Social Sciences, California 
% Institute of Technology, Pasadena, CA 91125.

% Rowe, D. B. and Press, S. J. (1998). Gibbs Sampling 
% and Hill Climbing in Bayesian Factor Analysis. Technical 
% Report No. 255 Department of Statistics, University of 
% California, Riverside.

% Press, S. J. and Shigemasu, K. (1997). Bayesian inference in 
% Factor Analysis-Revised. Technical Report No. 243, Department 
% of Statistics, University of California, Riverside.

% Press, S. J. and Shigemasu, K. (1989). Bayesian inference in 
% Factor Analysis. In Contributions to Probability and 
% Statistics, Chapter 15, Springer-Verlag.

% Press S. J. (1989). Bayesian Statistics: Principles Models, 
% and Applications, John Wiley and Sons, New York.

% Press, S. J. (1982). Applied Multivariate Analysis: Using 
% Bayesian and Frequentist Methods of Inference. Robert E. 
% Krieger Publishing Company, Malabar, Florida.

% Geman, S. and Geman, D. (1984). Stochastic relaxation,
% Gibbs distributions and the Bayesian restoration of images.
% IEEE Trans. on pattern analysis and machine intelligence
% 6, 721-741.

% Lindley, D. V. and Smith, A. F. M. (1972). Bayes estimates for 
% the linear model, Journal of the Royal Statistical Society B, 
% Volume 34, No. 1.

% O'Hagan, Anthony (1994). Kendall's Advanced Theory 
% of Statistics, Volume 2B Bayesian Inference. John Wiley 
% and Sons, Inc., New York.

06 07 02 05 08 07 08 08 03 08 09 07 05 07 10
09 10 05 08 10 09 09 10 05 09 09 08 08 08 10
07 08 03 06 09 08 09 07 04 09 09 08 06 08 10
05 06 08 05 06 05 09 02 08 04 05 08 07 06 05
06 08 08 08 04 05 09 02 08 05 05 08 08 07 07
07 07 07 06 08 07 10 05 09 06 05 08 06 06 06
09 09 08 08 08 08 08 08 10 08 10 08 09 08 10
09 09 09 08 09 09 08 08 10 09 10 09 09 09 10
09 09 07 08 08 08 08 05 09 08 09 08 08 08 10
04 07 10 02 10 10 07 10 03 10 10 10 09 03 10
04 07 10 00 10 08 03 09 05 09 10 08 10 02 05
04 07 10 04 10 10 07 08 02 08 08 10 10 03 07
06 09 08 10 05 04 09 04 04 04 05 04 07 06 08
08 09 08 09 06 03 08 02 05 02 06 06 07 05 06
04 08 08 07 05 04 10 02 07 05 03 06 06 04 06
06 09 06 07 08 09 08 09 08 08 07 06 08 06 10
08 07 07 07 09 05 08 06 06 07 08 06 06 07 08
06 08 08 04 08 08 06 04 03 03 06 07 02 06 04
06 07 08 04 07 08 05 04 04 02 06 08 03 05 04
04 08 07 08 08 09 10 05 02 06 07 09 08 08 09
03 08 06 08 08 08 10 05 03 06 07 08 08 05 08
09 08 07 08 09 10 10 10 03 10 08 10 08 10 08
07 10 07 09 09 09 10 10 03 09 09 10 09 10 08
09 08 07 10 08 10 10 10 02 09 07 09 09 10 08
06 09 07 07 04 05 09 03 02 04 04 04 04 05 04
07 08 07 08 05 04 08 02 03 04 05 06 05 05 06
02 10 07 09 08 09 10 05 03 05 06 07 06 04 05
06 03 05 03 05 03 05 00 00 03 03 00 00 05 00
04 03 04 03 03 00 00 00 00 04 04 00 00 05 00
04 06 05 06 09 04 10 03 01 03 03 02 02 07 03
05 05 04 07 08 04 10 03 02 05 05 03 04 08 03
03 03 05 07 07 09 10 03 02 05 03 07 05 05 02
02 03 05 07 07 09 10 03 02 02 03 06 04 05 02
03 04 06 04 03 03 08 01 01 03 03 03 02 05 02
06 07 04 03 03 00 09 00 01 00 02 03 01 05 03
09 08 05 05 06 06 08 02 02 02 04 05 06 06 03
04 09 06 04 10 08 08 09 01 03 09 07 05 03 02
04 09 06 06 09 09 07 09 01 02 10 08 05 05 02
10 06 09 10 09 10 10 10 10 10 08 10 10 10 10
10 06 09 10 09 10 10 10 10 10 10 10 10 10 10
10 07 08 00 02 01 02 00 10 02 00 03 00 00 10
10 03 08 00 01 01 00 00 10 00 00 00 00 00 10
03 04 09 08 02 04 05 03 06 02 01 03 03 03 08
07 07 07 06 09 08 08 06 08 08 10 08 08 06 05
09 06 10 09 07 07 10 02 01 05 05 07 08 04 05
09 08 10 10 07 09 10 03 01 05 07 09 09 04 04
00 07 10 03 05 00 10 00 00 02 02 00 00 00 00
00 06 10 01 05 00 10 00 00 02 02 00 00 00 00

% Copyright (c) 2002, Daniel B. Rowe
% This Matlab program performs the computations 
% for Bayesian Source Separation as described
% in Rowe 2002 with conjugate priors.
% There are references listed at the end.

clear
fignum=0;

% Have data arranged so that n is the number 
% of rows and p is the number of columns
% Each row is an observation vector.

load 'plankton.txt'   % Reads in data from plankton.txt
[nn,p]=size(plankton);
n0=nn/2;
n=nn/2;
for i=1:n0
X0(i,1:p)=plankton(2*(i-1)+1,1:p);
end
for i=1:n
X(i,1:p)=plankton(2*i,1:p);
end
clear plankton


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hyperparameter assessment                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Assessed Hyperparameters')

[pc,score,latent,tsquare] = princomp(X0);
S0=score(:,1:5)
m=5;

D=zeros(m+1,m+1);
D=inv([ones(n0,1),S0]'*[ones(n0,1),S0])
C0=zeros(p,m+1);
C0=X0'*[ones(n0,1),S0]*D
mu0=zeros(p,1);
Lambda0=zeros(p,m);
mu0=C0(:,1:1);
Lambda0=C0(:,2:m+1);

nu=n0               % nu must be a positive integer
Q=(X-[ones(n0,1),S0]*C0')'*(X-[ones(n0,1),S0]*C0') 

eta=n0;
V=S0'*S0;


% These numbers can be changed depending on the data set.
ICMiterations=25000;      % # of ICM iterations
Gibbsburn=10000;          % # of burn in variates
Gibbsiterations=25000;    % # of samples for estimates


% Do not change anything below here.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimation                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Dinv=zeros(m+1,m+1);
Dinv=inv(D);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ICM algorithm estimates                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize variables
SICM=zeros(n,m); 
ZICM=zeros(n,m+1);
CICM=zeros(p,m+1);
RICM=zeros(m,m);
SigmaICM=zeros(p,p);
ZICM=[ones(n,1),SICM];


% Initial value
SICM=randn(n,m);

for count=1:ICMiterations  
 ZICM=[ones(n,1),SICM];
 CICM=(X'*ZICM+C0*Dinv)*inv(Dinv+ZICM'*ZICM);
 SigmaICM=[(X-ZICM*CICM')'*(X-ZICM*CICM')+(CICM-C0)*Dinv*(CICM-C0)'+Q]/(n+nu+m+1);
 RICM=((SICM-S0)'*(SICM-S0)+V)/(n+eta);
 SICM1=S0*inv(RICM)+(X-ones(n,1)*CICM(:,1)')*inv(SigmaICM)*CICM(:,2:m+1); 
 SICM=SICM1*inv(inv(RICM)+CICM(:,2:m+1)'*inv(SigmaICM)*CICM(:,2:m+1)); 
end
 
 CICM                          
 SigmaICM
 RICM
 SICM

%%%ICM Significance
WICM=zeros(m+1,m+1);
WICM=inv(Dinv+ZICM'*ZICM);

disp(' ')
disp('Computing marginal conditional posterior statistics')
% Computing marginal t-statistics for each coefficient
zICM=zeros(p,m+1);
for col=1:m+1
   for row=1:p
    zICM(row,col)=CICM(row,col)/sqrt(WICM(col,col)*SigmaICM(row,row)) ; 
   end
end
zICM



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gibbs sampling algorithm estimates             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize variables 
% The starting values don't matter
% But convergence is faster with good values
SGibbs=SICM;
CGibbs=CICM;
RGibbs=RICM;
SigmaGibbs=SigmaICM;
ZGibbs=[ones(n,1),SGibbs];
WGibbs=zeros(m+1,m+1);
WGibbs=inv(Dinv+ZGibbs'*ZGibbs);


Stot=zeros(n,m);
Ctot=zeros(p,m+1);
C2tot=zeros(p*(m+1),p*(m+1));
Rtot=zeros(m,m);
Sigmatot=zeros(p,p);

YC = zeros(p,m+1);
YSig = zeros(n+nu+m+1+p+1,p);
SGibbs=zeros(n,m);
YR = zeros(n+eta+m+1,m);

for count=1:Gibbsburn+Gibbsiterations   
   
  MC=(C0*Dinv+X'*ZGibbs)*WGibbs;
  YC = randn(p,m+1);
  AtC=chol(SigmaGibbs); 
  BtC=chol(WGibbs);
  CGibbs=AtC'*YC*BtC+MC;
   
  YSig = randn(n+nu+m+1+p+1,p);
  AtSig=chol((X-ZGibbs*CGibbs')'*(X-ZGibbs*CGibbs')+(CGibbs-C0)*Dinv*(CGibbs-C0)'+Q);
  SigmaGibbs=AtSig'*inv(YSig'*YSig)*AtSig;
  
  YR = randn(n+eta+m+1,m);
  AtR=chol((SGibbs-S0)'*(SGibbs-S0)+V);
  RGibbs=AtR'*inv(YR'*YR)*AtR;
  
  MS1=S0*inv(RGibbs)+(X-ones(n,1)*CGibbs(:,1)')*inv(SigmaGibbs)*CGibbs(:,2:m+1);
  MS=MS1*inv(inv(RGibbs)+CGibbs(:,2:m+1)'*inv(SigmaGibbs)*CGibbs(:,2:m+1));
  SGibbs=randn(n,m);
  BtS=chol(inv(inv(RGibbs)+CGibbs(:,2:m+1)'*inv(SigmaGibbs)*CGibbs(:,2:m+1)));
  SGibbs=SGibbs*BtS+MS;
  
  ZGibbs=[ones(n,1),SGibbs];
  WGibbs=inv(Dinv+ZGibbs'*ZGibbs);
  
  if count>Gibbsburn
   Stot=Stot+SGibbs;
   Ctot=Ctot+CGibbs;
   C2tot=C2tot+vec(CGibbs)*vec(CGibbs)';
   Rtot=Rtot+RGibbs;
   Sigmatot=Sigmatot+SigmaGibbs;
  end
end


SGibbs=Stot/Gibbsiterations
CGibbs=Ctot/Gibbsiterations
DeltaGibbs=C2tot/Gibbsiterations-vec(CGibbs)*vec(CGibbs)';
RGibbs=Rtot/Gibbsiterations
SigmaGibbs=Sigmatot/Gibbsiterations 

%%%Gibbs Significance
disp(' ')
disp('Computing marginal posterior statistics')
% Computing marginal z-statistics for each coefficient
zGibbs=zeros(p,m+1);
for col=1:m+1
   for row=1:p
    zGibbs(row,col)=CGibbs(row,col)/sqrt(DeltaGibbs(p*(col-1)+row,p*(col-1)+row)); 
   end
end
zGibbs

fignum=fignum+1;
figure(fignum)
 subplot(9,2,1)
  plot(1:n,X(:,1),'k')
  axis tight
 subplot(9,2,3)
  plot(1:n,X(:,2),'k')
  axis tight
 subplot(9,2,5)
  plot(1:n,X(:,3),'k')
  axis tight
 subplot(9,2,7)
  plot(1:n,X(:,4),'k')
  axis tight
 subplot(9,2,9)
  plot(1:n,X(:,5),'k')
  axis tight
 subplot(9,2,11)
  plot(1:n,X(:,6),'k')
  axis tight
 subplot(9,2,13)
  plot(1:n,X(:,7),'k')
  axis tight
 subplot(9,2,15)
  plot(1:n,X(:,8),'k')
  axis tight
 subplot(9,2,17)
 plot(1:n,X(:,9),'k')
 axis tight
 
 subplot(9,2,6)
  plot(1:n,SGibbs(:,1),'k--',1:n,SICM(:,1),'k:')
  axis tight
 subplot(9,2,8)
  plot(1:n,SGibbs(:,2),'k--',1:n,SICM(:,2),'k:')
  axis tight
 subplot(9,2,10)
  plot(1:n,SGibbs(:,3),'k--',1:n,SICM(:,3),'k:')
  axis tight
 subplot(9,2,12)
  plot(1:n,SGibbs(:,4),'k--',1:n,SICM(:,4),'k:')
  axis tight
 subplot(9,2,14)
  plot(1:n,SGibbs(:,5),'k--',1:n,SICM(:,5),'k:')
  axis tight
  

%save bssnat


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% References                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Rowe, D.B. (2003). Multivariate Bayesian Statistics: 
% Models for Source Separation and Signal Unmixing.
% CRC Press, Boca Raton, Florida.

% Rowe, D.B. (2003). Bayesian Source Separation of 
% Functional Sources. Journal of Interdisciplinary 
% Mathematics, Vol. 6, No. 2, pp. 129-138.

% Rowe, D.B. (2002). Jointly Distributed Mean and 
% Mixing Coefficients for Bayesian Source Separation 
% using MCMC and ICM. Monte Carlo Methods and Applications, 
% Vol. 8, No. 4, pp. 395-403.

% Rowe, D.B. (2001). Bayesian Source Separation for 
% Reference Function Determination in fMRI. Magnetic 
% Resonance in Medicine, Volume 46, Issue 2, pp. 374-378.

% Rowe, D.B. (2002). A Bayesian Approach to Blind 
% Source Separation. Journal of Interdisciplinary 
% Mathematics, Vol. 5, No. 1, pp. 49-76. 

% Rowe, D.B. (2001). Bayesian Source Separation with 
% Jointly Distributed Mean and Mixing Coefficients 
% via MCMC and ICM. Social Science Working Paper 1119, 
% Division of Humanities and Social Sciences, California 
% Institute of Technology, Pasadena, CA 91125. 

% Rowe, D.B. (2001). A Model for Bayesian Source 
% Separation With The Overall Mean. Social Science 
% Working Paper 1118, Division of Humanities and Social 
% Sciences, California Institute of Technology, Pasadena, 
% CA 91125. 

% Rowe, D. B. and Press, S. J. (1998). Gibbs Sampling 
% and Hill Climbing in Bayesian Factor Analysis. 
% Technical Report No. 255 Department of Statistics, 
% University of California, Riverside.

% Geman, S. and Geman, D. (1984). Stochastic relaxation,
% Gibbs distributions and the Bayesian restoration of images.
% IEEE Trans. on pattern analysis and machine intelligence
% 6, 721-741.

% Lindley, D. V. and Smith, A. F. M. (1972). Bayes estimates  
% for the linear model, Journal of the Royal Statistical  
% Society B, Volume 34, No. 1.

% O'Hagan, Anthony (1994). Kendall's Advanced Theory 
% of Statistics, Volume 2B Bayesian Inference. John Wiley 
% and Sons, Inc., New York.

 1.792  0.489 43.485  0.814 25.570  0.651  0.163  0.000  0.163
 3.203  0.712 37.722  0.356 30.961  0.712  0.356  0.000  0.000
 2.364  1.709 47.009  0.855 20.513  1.709  1.282  0.427  0.000
 1.124  0.562 47.191  1.124 12.360  2.247  3.933  0.562  0.562
 0.671  1.007 43.624  3.020 15.436  1.007  0.336  0.671  0.336
 1.149  0.766 52.874  0.766 12.261  0.000  0.383  2.299  0.000
 1.990  0.498 53.234  3.980  6.965  0.000  0.498  0.995  0.000
 2.222  2.222 45.926  2.222 13.333  2.963  1.481  1.481  1.481
 1.786  1.190 49.405  1.786 10.714  1.786  0.595  0.595  0.000
 0.621  0.621 36.025  2.484 10.519  0.621  1.242  1.863  0.000
 1.418  0.000 46.099  2.837  9.220  4.255  0.709  2.836  0.000
 0.000  0.000 38.298  0.709 11.348  2.837  1.418  5.674  0.000
 0.498  0.498 48.756  0.000  5.970  1.990  0.498  2.985  0.000
 1.379  1.034 42.069  0.690  8.621  2.069  2.759  1.724  0.690
 0.662  0.000 46.358  0.000 11.921  0.000  1.987  3.311  0.000
 3.429  1.143 45.714  1.143 14.286  1.714  0.571  3.429  0.571
 2.899  2.899 42.995  0.000 14.010  1.449  2.415  2.415  0.483
 1.198  1.796 50.299  1.198  8.383  2.994  0.599  0.599  0.599
 1.887  2.516 38.994  3.145  7.547  2.516  1.258  1.258  0.000
 5.143  2.857 38.286  0.000 13.714  1.143  1.143  1.143  0.000
 3.067  0.613 37.423  1.227 13.497  2.761  1.227  0.000  0.307
 1.961  2.614 41.830  3.268 11.765  1.307  1.307  0.654  0.000
 1.515  2.020 37.374  1.010 12.626  2.020  0.000  0.505  0.000
 1.422  2.844 38.389  1.422 16.114  0.948  0.000  0.474  0.000
 1.630  1.630 36.957  2.174 10.870  2.174  0.000  0.000  0.000
 1.571  1.571 37.696  1.571 10.995  4.188  2.094  2.618  1.047
 1.826  3.196 36.073  0.913 12.329  2.283  0.457  0.913  0.457
 0.926  3.241 28.241  0.463 12.037  0.926  0.463  1.852  0.463
 1.379  2.414 35.517  0.345 11.679  0.345  0.000  4.828  0.000
 1.036  6.218 34.197  1.036 14.508  0.518  0.000  1.554  0.518
 0.649  3.896 39.610  3.896 13.636  1.299  0.543  0.649  0.000
 1.485  7.426 29.208  2.475 15.842  1.485  2.970  1.485  0.000
 1.087  0.000 42.391  1.630 15.761  1.630  2.174  1.087  0.000
 3.404  0.426 32.766  4.255 13.191  2.128  3.830  0.851  1.700
 1.429  0.476 42.381  2.857 10.952  1.905  0.476  0.952  1.900
 1.449  3.623 36.957  0.000 15.942  3.623  0.725  1.449  0.720
 1.685  1.685 48.315  2.809 10.674  1.124  1.124  1.124  0.000
 0.772  0.386 40.927  0.772 15.444  2.703  0.000  0.772  0.380
 1.266  1.266 37.975  2.532 18.143  3.376  2.110  0.422  0.000
 3.627  0.518 41.451  1.554 16.580  0.518  2.591  1.554  0.000
 1.869  1.402 37.850  2.804 12.617  2.336  9.813  0.467  0.930
 3.509  2.456 42.105  2.105 12.281  1.053  2.456  0.000  0.000
 0.904  0.904 44.578  1.205 14.759  0.602  1.506  0.602  0.000
 1.449  0.483 43.961  3.865 12.560  1.449  2.899  0.000  0.000
 1.299  0.649 38.961  0.325 17.208  1.945  4.545  1.948  0.000
 0.000  0.741 33.333  2.222 22.222  2.222  0.741  0.000  0.000
 2.513  4.523 35.176  1.005 20.603  0.000  0.000  0.000  0.000
 1.026  0.513 42.051  2.051 16.410  2.051  0.513  2.051  0.000
 0.565  0.565 44.068  3.955 10.169  1.695  9.605  3.390  0.000
 1.523  0.000 34.518  2.030 20.305  2.030  1.523  1.015  0.000
 0.508  0.000 40.609  0.508 21.827  0.508  3.046  0.000  0.000
 0.000  2.703 28.649  1.622 24.324  3.784  2.162  3.243  0.000
 0.629  4.403 39.623  0.629 10.063  3.145  5.660  5.031  0.000
 0.800  2.400 50.400  1.600 11.200  2.400  4.800  0.000  0.000
 1.630  0.543 54.348  2.174  7.609  3.804  1.630  2.717  0.000
 0.000  0.543 32.609  1.087 11.413  4.891  3.804  2.717  0.000
 1.622  1.081 32.973  2.162 11.892  3.784  9.780  0.541  0.000
 1.762  0.000 33.921  0.000 16.740  2.643  9.251  2.643  0.000
 1.418  0.000 36.879  0.709 11.348  4.255  4.965  4.965  0.709
 1.136  2.841 49.432  2.273 11.932  2.273  0.568  0.000  0.000
 0.893  3.561 33.036  5.357 13.393  2.679  4.464  0.893  0.893
 3.636  1.212 35.758  2.424  6.061  6.061  3.030  0.000  0.000
 3.448  1.478 29.064  3.448 14.778  4.433  2.955  0.000  0.000
 1.342  2.685 34.228  3.356 12.081  2.685  2.685  4.027  0.000
 4.435  2.419 33.468  0.806 17.742  3.226  0.000  4.032  0.000
 2.158  2.158 34.532  2.158 15.826  5.036  0.719  2.158  0.000
 0.000  4.545 38.636  0.000 15.152  1.515  2.273  2.273  0.758
 1.235  0.000 41.975  0.000 12.346  1.852  0.617  2.469  0.000
 1.508  1.508 38.191  0.503  3.518  1.508  1.508  2.010  0.503
 3.550  2.367 47.337  2.367  5.917 10.059  0.000  0.592  0.000
 5.344  0.000 39.695  1.527 13.740  6.870  0.763  0.000  0.000
 5.455  0.606 43.636  1.818 10.303  7.273  0.605  0.000  0.000
 0.000  0.000 38.095  3.571  4.762  9.524  3.571  0.000  1.190
 2.609  1.304 33.043  1.739  9.130  3.913  3.478  0.435  0.000
 1.604  1.604 33.690  0.000 19.251  2.139  3.209  3.209  0.535
 1.899  0.000 34.177  2.532 12.025  4.430  2.532  1.266  0.000
 2.041  0.816 36.327  2.041 20.000  2.449  2.449  1.224  0.408
 0.595  2.976 50.000  0.000  7.738  6.548  2.381  0.595  0.000
 0.000  6.130 35.249  0.000 10.728  0.000  0.383  0.383  0.000
 0.372  5.576 37.918  0.372 15.613  0.743  0.000  0.372  0.000
 3.582  5.373 38.209  0.896 17.015  0.896  0.000  0.896  0.299
 2.362  2.362 36.220  3.150 14.173  1.969  0.787  1.575  0.000
 2.105  4.211 26.842  1.053 13.684  4.737  5.263  2.105  0.000
 2.381  3.175 32.143  1.190 17.460  1.587  0.397  1.190  0.000
 0.455  0.909 37.273  0.455 24.091  3.182  0.455  0.455  0.909
 0.858  3.863 31.760  1.717 21.888  7.296  4.721  0.858  0.000
 2.769  1.231 43.385  1.231  2.769  4.000  6.462  3.077  0.000
 0.658  1.316 52.632  0.000  3.289  1.974  3.947  0.658  0.000
 3.448  0.575 35.632  1.149 14.368  0.000  4.598  0.575  0.000
 1.689  0.676 26.689  2.027  8.108  4.392 13.176  2.027  1.689
 1.533  0.000 35.249  0.383  9.195  2.682 13.793  1.533  0.000
 1.064  0.000 40.957  1.596  6.915  2.660  3.723  2.660  0.000
 1.394  0.348 36.585  1.045  8.014  3.833  6.969  1.394  0.000
 0.000  0.000 35.533  1.015 13.706  7.614  3.553  0.000  0.000
 1.970  2.463 39.901  0.493 15.764  3.941  0.985  0.493  0.493
 1.471  2.206 34.559  2.941 15.441  1.471  0.000  0.735  0.000
 1.613  0.403 42.742  1.210 16.129  2.823  2.823  0.403  0.000
 0.000  0.498 44.776  2.488 19.900  0.995  1.990  0.995  0.498
 0.448  0.448 40.359  4.484 12.556  2.242  6.278  0.897  0.000
 2.717  0.000 32.065  3.261 15.761  1.087  6.522  1.087  0.000
 1.887  1.887 34.906  1.415 12.264  1.415  3.302  1.415  0.472
 1.342  2.013 24.161  3.356 11.409  1.342  9.396  0.000  0.671
 1.633  0.816 24.898  2.449  6.531  0.408 12.245  2.041  0.000
 1.548  0.310 31.269  1.548  9.288  0.000  9.288  4.644  0.000
 1.093  0.546 31.694  1.639 14.208  0.000 19.672  4.372  0.000
 2.183  1.747 33.188  0.437 13.974  0.437  4.367  1.747  1.747
 1.878  0.469 24.883  1.878 14.085  1.408  9.390  0.939  0.000
 2.286  2.286 37.143  1.714  8.000  1.714  8.000  4.571  0.000
 3.911  2.793 32.961  1.117 14.525  1.117  2.793  0.559  0.000
 0.658  0.658 34.868  4.605 15.789  1.316  3.947  1.974  0.000

% Copyright (c) 2002, Daniel B. Rowe
% This Matlab program performs the computations 
% for Bayesian Source Separation with unobservable
% and observable sources as described
% in Rowe 2002 with conjugate priors.
% There are references listed at the end.

clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Reads in data from fmrisim.txt               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Have data arranged so that n is the number 
% of rows and p is the number of columns
% Each row is an observation vector.

n=128;
p=16;
q=1;
fignum=0;

% Reads in data from bssuodata.txt 
disp(' ')
disp('Reading in the data')
load 'fmrisim.txt'

disp(' ')
disp('The data matrix X is') 
X=fmrisim(:,1:p)
X0=X;      % Empirical Bayesian Approach

disp(' ')
disp('The design matrix U is')
U=[ones(n,1),(1:n)']      %could be read from a file

clear fmrisim



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Hyperparameter assessment                      %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Assessed Hyperparameters')

m=2;
S0=zeros(n,m); 
S0temp=[1,    -1;
        1,    -1;
        1,    -1;
        1,    -1;
        1,    -1;
        0,     0;
       -1,     1;
       -1,     1;
       -1,    -1;
       -1,    -1;
       -1,    -1;
       -1,    -1;
       -1,    -1;
       -1,    -1;
       -1,    -1;
       -1,    -1];
S0=[S0temp;S0temp;S0temp;S0temp;S0temp;S0temp;S0temp;S0temp];
clear S0temp


D=zeros(m+q+1,m+q+1);
D=inv([U,S0]'*[U,S0])
C0=zeros(p,m+q+1);
C0=X0'*[U,S0]*D
B0=zeros(p,q+1);
Lambda0=zeros(p,m);
B0=C0(:,1:q+1);
Lambda0=C0(:,q+2:m+q+1);

nu=n               % nu must be a positive integer
Q=(X0-[U,S0]*C0')'*(X0-[U,S0]*C0')


%%% Regression Significance i.e. t-statistics
zcrit=5;
tcrit=zcrit;

treg1=zeros(p,1);
treg2=zeros(p,1);

for j=1:p
   treg1(j,1)=C0(j,q+1+1)/sqrt(Q(j,j)/(n-m-q-1)*D(q+1+1,q+1+1));
   treg2(j,1)=C0(j,q+1+2)/sqrt(Q(j,j)/(n-m-q-1)*D(q+1+2,q+1+2));
end
treg1
treg2

xdim=4;
ydim=4;
actreg1=reshape(treg1,xdim,ydim)
actreg2=reshape(treg2,xdim,ydim)

mycmap=jet(16);        % Changing colormap
mymap(1,:)=mycmap(7,:);
mymap(2,:)=mycmap(6,:);
mymap(3,:)=mycmap(5,:);
mymap(4,:)=mycmap(4,:);
mymap(5,:)=mycmap(3,:);
mymap(6,:)=mycmap(2,:);
mymap(7,:)=mycmap(1,:);
mymap(10,:)=mycmap(16,:);
mymap(11,:)=mycmap(15,:);
mymap(12,:)=mycmap(14,:);
mymap(13,:)=mycmap(13,:);
mymap(14,:)=mycmap(12,:);
mymap(15,:)=mycmap(11,:);
mymap(16,:)=mycmap(10,:);
mymap(8:9,:)=1;
mymap;

fignum=fignum+1;
figure(fignum)
colormap(mymap)
imagesc(actreg1,[-15,15])
axis equal
axis tight
caxis([-15,15]);
title('Unthresholded Reg 1')
colorbar

actreg1t=zeros(xdim,ydim);

% Thresholding
for count1=1:xdim,
   for count2=1:ydim,
      if (actreg1(count1,count2)) > tcrit
         actreg1t(count1,count2)=actreg1(count1,count2);
      elseif (actreg1(count1,count2)) < -tcrit
         actreg1t(count1,count2)=actreg1(count1,count2);
      else
         actreg1t(count1,count2)=0.0;
      end
   end
end
actreg1t

fignum=fignum+1;
figure(fignum)
colormap(mymap)
imagesc(actreg1t,[-15,15])
axis equal
axis tight
title('Thresholded Reg 1')
colorbar

actreg2t=zeros(xdim,ydim);
fignum=fignum+1;
figure(fignum)
colormap(mymap)
imagesc(actreg2,[-15,15])
axis equal
axis tight
caxis([-15,15]);
title('Unthresholded Reg 2')
colorbar

% Thresholding
for count1=1:xdim,
   for count2=1:ydim,
      if actreg2(count1,count2)>tcrit
         actreg2t(count1,count2)=actreg2(count1,count2);
      elseif (actreg2(count1,count2)) < -tcrit
         actreg2t(count1,count2)=actreg2(count1,count2);
      else
         actreg2t(count1,count2)=0.0;
      end
   end
end
actreg2t

fignum=fignum+1;
figure(fignum)
colormap(mymap)
imagesc(actreg2t,[-15,15])
axis equal
axis tight
title('Thresholded Reg 2')
colorbar

%%%


er=.0625;
varr=.25;          %sqrt(.0625)=.25
eta=er*er/varr/2+6;
eta=eta-rem(eta,1)
v0=er*(eta-4);
V=v0*eye(m)




% These numbers can be changed depending on the data set.
ICMiterations=25000;      % # of ICM iterations
Gibbsburn=5000;          % # of burn in samples
Gibbsiterations=25000;   % # of samples for estimates


% Do not change any thing below here.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimation                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Dinv=zeros(m+q+1,m+q+1);
Dinv=inv(D);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ICM algorithm estimates                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize variables
SICM=zeros(n,m); 
ZICM=zeros(n,m+q+1);
CICM=zeros(p,m+q+1);
RICM=zeros(m,m);
SigmaICM=zeros(p,p);

% Initial value
SICM=S0;
RICM=V/eta;

for count=1:ICMiterations  
 ZICM=[U,SICM];
 CICM=(X'*ZICM+C0*Dinv)*inv(Dinv+ZICM'*ZICM);
 SigmaICM=[(X-ZICM*CICM')'*(X-ZICM*CICM')+(CICM-C0)*Dinv*(CICM-C0)'+Q]/(n+nu+m+1);
 SICM1=S0*inv(RICM)+(X-U*CICM(:,1:q+1)')*inv(SigmaICM)*CICM(:,q+2:m+q+1);      %
 
 SICM=SICM1*inv(inv(RICM)+CICM(:,q+2:m+q+1)'*inv(SigmaICM)*CICM(:,q+2:m+q+1));
 RICM=((SICM-S0)'*(SICM-S0)+V)/(n+eta);
end
 
 CICM                          
 SigmaICM
 RICM
 SICMclear
 

%%%ICM Significance
W=inv(Dinv+[U,SICM]'*[U,SICM]);

zICM1=zeros(p,1);
zICM2=zeros(p,1);
for j=1:p
   zICM1(j,1)=CICM(j,q+1+1)/sqrt(W(q+1+1,q+1+1)*SigmaICM(j,j));
   zICM2(j,1)=CICM(j,q+1+2)/sqrt(W(q+1+2,q+1+2)*SigmaICM(j,j));
end
zICM1
zICM2

actICM1=reshape(zICM1,xdim,ydim)
actICM2=reshape(zICM2,xdim,ydim)

fignum=fignum+1;
figure(fignum)
colormap(mymap)
imagesc(actICM1,[-15,15])
axis equal
axis tight
caxis([-15,15]);
title('Unthresholded ICM 1')
colorbar

actICM1t=zeros(xdim,ydim);

% Thresholding
for count1=1:xdim,
   for count2=1:ydim,
      if (actICM1(count1,count2)) > zcrit
         actICM1t(count1,count2)=actICM1(count1,count2);
      elseif (actICM1(count1,count2)) < -zcrit
         actICM1t(count1,count2)=actICM1(count1,count2);
      else
         actICM1t(count1,count2)=0.0;
      end
   end
end
actICM1t

fignum=fignum+1;
figure(fignum)
colormap(mymap)
imagesc(actICM1t,[-15,15])
axis equal
axis tight
caxis([-15,15]);
title('Thresholded ICM 1')
colorbar

actICM2t=zeros(xdim,ydim);
fignum=fignum+1;
figure(fignum)
colormap(mymap)
imagesc(actICM2,[-15,15])
axis equal
axis tight
caxis([-15,15]);
title('Unthresholded ICM 2')
colorbar

% Thresholding
for count1=1:xdim,
   for count2=1:ydim,
      if actICM2(count1,count2)>zcrit
         actICM2t(count1,count2)=actICM2(count1,count2);
      elseif (actICM2(count1,count2)) < -zcrit
         actICM2t(count1,count2)=actICM2(count1,count2);
      else
         actICM2t(count1,count2)=0.0;
      end
   end
end
actICM2t


fignum=fignum+1;
figure(fignum)
colormap(mymap)
imagesc(actICM2t,[-15,15])
axis equal
axis tight
title('Thresholded ICM 2')
colorbar

%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gibbs sampling algorithm estimates             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initialize variables 
% The starting values don't matter
% But convergence is faster with good values
SGibbs=SICM;
CGibbs=CICM;
RGibbs=RICM;
SigmaGibbs=SigmaICM;
ZGibbs=[U,SGibbs];
WGibbs=zeros(m+q+1,m+q+1);
WGibbs=inv(Dinv+ZGibbs'*ZGibbs);
SigmaGibbsinv=zeros(p,p);
RGibbsinv=RGibbs(m,m);

Stot=zeros(n,m);
Ctot=zeros(p,m+q+1);
C2tot=zeros(p*(m+q+1),p*(m+q+1));
Rtot=zeros(m,m);
Sigmatot=zeros(p,p);

YC = zeros(p,m+q+1);
YP = zeros(n+nu+m+1+p+1,p);
SGibbs=zeros(n,m);
YR= zeros(n+eta+m+1,m);


for count=1:Gibbsburn+Gibbsiterations   
   
  MC=(C0*Dinv+X'*ZGibbs)*WGibbs;
  YC = randn(p,m+q+1);
  AtC=chol(SigmaGibbs); 
  BtC=chol(WGibbs);
  CGibbs=AtC'*YC*BtC+MC;
   
  YSig = randn(n+nu+m+1+p+1,p);
  AtSig=chol((X-ZGibbs*CGibbs')'*(X-ZGibbs*CGibbs')+(CGibbs-C0)*Dinv*(CGibbs-C0)'+Q);
  SigmaGibbs=AtSig'*inv(YSig'*YSig)*AtSig;
  SigmaGibbsinv=inv(SigmaGibbs);
  
  YR = randn(n+eta+m+1,m);
  AR=chol((SGibbs-S0)'*(SGibbs-S0)+V);
  RGibbs=AR'*inv(YR'*YR)*AR;
  RGibbsinv=inv(RGibbs);
  
  MS1=S0*RGibbsinv+(X-U*CGibbs(:,1:q+1)')*SigmaGibbsinv*CGibbs(:,q+2:m+q+1);
  MS=MS1*inv(RGibbsinv+CGibbs(:,q+2:m+q+1)'*SigmaGibbsinv*CGibbs(:,q+2:m+q+1));
  SGibbs=randn(n,m);
  BtS=chol(inv(RGibbsinv+CGibbs(:,q+2:m+q+1)'*SigmaGibbsinv*CGibbs(:,q+2:m+q+1)));
  SGibbs=SGibbs*BtS+MS;
  
  ZGibbs=[U,SGibbs];
  WGibbs=inv(Dinv+ZGibbs'*ZGibbs);
  
  if count>Gibbsburn
   Stot=Stot+SGibbs;
   Ctot=Ctot+CGibbs;
   C2tot=C2tot+vec(CGibbs)*vec(CGibbs)';
   Rtot=Rtot+RGibbs;
   Sigmatot=Sigmatot+SigmaGibbs;
  end
  
end


SGibbs=Stot/Gibbsiterations
CGibbs=Ctot/Gibbsiterations
DeltaGibbs=C2tot/Gibbsiterations-vec(CGibbs)*vec(CGibbs)';
RGibbs=Rtot/Gibbsiterations
SigmaGibbs=Sigmatot/Gibbsiterations 


%%%Gibbs Significance
UpsilonGibbs=DeltaGibbs(p*(q+1)+1:p*(q+1+m),p*(q+1)+1:p*(q+1+m));

zGibbs1=zeros(p,1);
zGibbs2=zeros(p,1);
for j=1:p
   zGibbs1(j,1)=CGibbs(j,q+1+1)/sqrt(UpsilonGibbs(j,j));
   zGibbs2(j,1)=CGibbs(j,q+1+2)/sqrt(UpsilonGibbs(p+j,p+j));
end
zGibbs1
zGibbs2

actGibbs1=reshape(zGibbs1,xdim,ydim)
actGibbs2=reshape(zGibbs2,xdim,ydim)

fignum=fignum+1;
figure(fignum)
colormap(mymap)
imagesc(actGibbs1,[-15,15])
axis equal
axis tight
caxis([-15,15]);
title('Unthresholded Gibbs 1')
colorbar

actGibbs1t=zeros(xdim,ydim);

% Thresholding
for count1=1:xdim,
   for count2=1:ydim,
      if (actGibbs1(count1,count2)) > zcrit
         actGibbs1t(count1,count2)=actGibbs1(count1,count2);
      elseif (actICM1(count1,count2)) < -zcrit
         actGibbs1t(count1,count2)=actGibbs1(count1,count2);
      else
         actGibbs1t(count1,count2)=0.0;
      end
   end
end
actGibbs1t

fignum=fignum+1;
figure(fignum)
colormap(mymap)
imagesc(actGibbs1t,[-15,15])
axis equal
axis tight
title('Thresholded Gibbs 1')
colorbar

actGibbs2t=zeros(xdim,ydim);
fignum=fignum+1;
figure(fignum)
colormap(mymap)
imagesc(actGibbs2,[-15,15])
axis equal
axis tight
caxis([-15,15]);
title('Unthresholded Gibbs 2')
colorbar

% Thresholding
for count1=1:xdim,
   for count2=1:ydim,
      if actGibbs2(count1,count2)>zcrit
         actGibbs2t(count1,count2)=actGibbs2(count1,count2);
      elseif (actGibbs2(count1,count2)) < -zcrit
         actGibbs2t(count1,count2)=actGibbs2(count1,count2);
      else
         actGibbs2t(count1,count2)=0.0;
      end
   end
end
actGibbs2t

fignum=fignum+1;
figure(fignum)
colormap(mymap)
imagesc(actGibbs2t,[-15,15])
axis equal
axis tight
title('Thresholded Gibbs 2')
colorbar
%%%



fignum=fignum+1;
figure(fignum)
plot((1:n),S0(1:n,1),'k-',1:n,SICM(:,1),'b:',1:n,SGibbs(:,1),'r--')
axis tight
fignum=fignum+1;
figure(fignum)
plot((1:n),S0(1:n,2),'k-',1:n,SICM(:,2),'b:',1:n,SGibbs(:,2),'r--')
axis tight







%save bssuonat


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% References                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Rowe, D.B. (2003). Multivariate Bayesian Statistics: 
% Models for Source Separation and Signal Unmixing.
% CRC Press, Boca Raton, Florida.

% Rowe, D.B. (2003). Bayesian Source Separation of 
% Functional Sources. Journal of Interdisciplinary 
% Mathematics, Vol. 6, No. 2, pp. 129-138.

% Rowe, D.B. (2002). Jointly Distributed Mean and 
% Mixing Coefficients for Bayesian Source Separation 
% using MCMC and ICM. Monte Carlo Methods and Applications, 
% Vol. 8, No. 4, pp. 395-403.

% Rowe, D.B. (2002). A Bayesian Approach to Blind 
% Source Separation. Journal of Interdisciplinary 
% Mathematics, Vol. 5, No. 1, 49-76. 

% Rowe, D.B. (2001). A Bayesian Unobservable/Observable 
% Source Separation Model and Activation Determination in 
% FMRI. Social Science Working Paper 1120, Division of 
% Humanities and Social Sciences, California Institute 
% of Technology, Pasadena, CA 91125.

% Rowe, D.B. (2001). Bayesian Source Separation for 
% Reference Function Determination in fMRI. Magnetic 
% Resonance in Medicine, Volume 46, Issue 2, 374-378.

% Geman, S. and Geman, D. (1984). Stochastic relaxation,
% Gibbs distributions and the Bayesian restoration of images.
% IEEE Trans. on pattern analysis and machine intelligence
% 6, 721-741.

% Lindley, D. V. and Smith, A. F. M. (1972). Bayes estimates for 
% the linear model, Journal of the Royal Statistical Society B, 
% Volume 34, No. 1.

% O'Hagan, Anthony (1994). Kendall's Advanced Theory 
% of Statistics, Volume 2B Bayesian Inference. John Wiley 
% and Sons, Inc., New York.

    4.9010    1.9994  -21.9189  -21.7303   -4.0009   -3.7932    7.0807  -15.6352  -15.1510   -4.0514   -3.7105   -7.8278   -8.3415  -12.9287   -6.2728    6.0848
   -4.5734   -3.9538    3.6736  -13.1921    3.4987    5.4460   -5.4218  -15.6906   -8.8049    8.6579  -13.3278   -8.4276  -18.0661  -11.4330  -13.4044   19.0082
   32.4915    4.4512  -14.3858  -33.7289    0.9312   20.5609    6.7420   -8.1803  -15.2093    2.5940   23.4415    5.1717   -5.1756  -14.7450   -0.7948   11.8997
   14.7543   15.9958  -24.4645  -30.7529    8.1685    8.6984    3.2014   -8.3559   -6.9913  -11.4789   11.2435    9.7663   -9.5264   10.2181  -21.6842   23.3477
   20.7249    9.7098  -19.5218  -20.2393   -0.4193   10.5938   -9.5399  -24.5124  -24.1702   -5.0983   15.0896   18.0689  -14.5524    3.5058    3.9580   -0.9261
    8.4978   -4.7917    4.5738   10.3824   -1.7374    6.5217   -0.2515   16.2873    0.4122   -0.1934    0.7659   27.2197   -1.5904   12.9799  -10.9977   16.0102
  -13.6645    2.0840   27.1155   30.8023  -14.6891   15.5099    4.7836   27.5905   18.6970    7.6175   -6.6469    2.9658   10.7529   20.9284    6.4699   -1.3337
  -23.7919   -0.1470   27.8385   23.5684    9.9169  -19.8999    8.0344   -1.1422    5.6786   -3.2137  -21.4698   -1.4841   20.7266   39.0278   -2.3342  -26.3885
   -3.4825    3.8653   22.8534   21.2419   -3.0367   -8.6905   -2.3251   11.3237    1.4902    1.8506  -19.7233    7.4224   20.2485   17.3196    3.6222  -23.1399
  -13.0906   -5.2907    0.6867  -13.2003    2.8526   -6.6459   13.4113   -5.7390   -6.2853    6.6102  -29.5579    5.6802    6.3349  -25.9328    2.3281  -19.6655
  -22.3610  -21.2362   -8.9510  -26.8911   -9.0925    9.4646  -12.7434    0.1410    0.8946   -0.9758   -9.8397    3.5593   -2.0347    3.2907   -8.7016   20.7414
  -20.4927   -8.3859   -4.5011  -35.3972    7.6256  -27.5008   28.3517   -1.0928  -26.3799   23.3284  -18.7032   -0.5509    1.7462  -38.6898    2.0116  -16.2805
   -4.2005   -4.5833  -26.9048  -10.9460   20.8466  -15.4052   -1.6710   -4.3890  -15.2324   17.3162    2.2344    9.9598  -23.0031   10.7452    0.3397   -7.6619
  -17.9032   11.1334  -11.5001  -20.5781   -2.3802   -3.5703    9.4844   -4.4165  -11.9579   -4.7436  -28.6248   14.2727  -12.2579  -44.7247   -8.4672  -15.0735
   -5.8635    9.7893  -16.9753  -14.3694    9.0918    4.6537    2.0451  -16.6763   -8.8167   -9.4683  -11.8010    3.7463  -25.0767  -18.3031   -7.3279    5.8507
   -3.0197   11.0039   -1.5644   10.6664   -2.6813   -9.9589   -1.6028   -3.3265    9.7403   10.9216   10.6828    6.8799   12.2380   -3.5256   17.2747   14.9835
    9.3862    1.4375   -6.0317  -20.6896   19.0875   15.6301   12.8910   -4.4635   -7.0832   13.2665    7.5973    4.6349  -15.8338   -9.9870   17.1211    9.7659
   11.2091  -16.7068   17.4870   -5.7415   13.8106   12.9594    5.6305    2.5079  -17.2860   -1.0994   17.6456   -5.5201  -22.1707    6.6872    3.7989   22.1766
   24.0692   16.3161    4.1110    0.7268   23.4883   16.5968    9.5680    2.7085   -6.3399    2.7141   32.7799   12.4842  -25.0780  -19.3914    2.1410   34.9045
   37.0524   -4.5770   -0.7774  -15.7181    7.6833   43.7452  -14.0794    1.0960  -12.5400    5.5197   33.5544   18.3134  -16.3090  -11.5120    5.9249   40.5421
   26.2432    5.1789   23.5658   -3.3308   19.6689   43.2871    5.5120   -4.4646   -1.5061    4.1549   17.9784   39.8036  -15.7159  -11.7877   22.9222   18.1655
    7.5548   -9.8672   21.9363    5.3114    0.6356    9.6248    7.8965   15.1175   -1.3432   20.0211   11.2008    3.3465    0.9853    3.4510    2.7360   19.1172
    7.6616   -4.2552   34.5340   16.5889   24.8097   -0.2729   18.5141   39.5829   23.5307   16.2689   -1.9567   -9.3266   15.5348   -2.1233   12.6464   19.2179
  -13.7756  -14.6162   29.8539   20.3289   16.3187    0.5135   11.1723   21.1792   14.2733    4.9921    5.8501   10.1894   20.9542   34.9027   -0.2531    2.1238
   25.6713   -8.2469   25.6658   22.5996    1.4852    1.7549   10.0044   26.0218   25.4294    7.6457   -4.4013   34.6161   19.0716   12.2552   17.7875   -6.1204
   -2.8557   17.2460    6.1349    8.0497   29.0066   14.9539   14.5111    7.5960   -3.1381   10.4805   14.3842   -2.2724   -8.9396   -4.6515    2.8486   15.0115
   -2.1552   13.4586   25.1275   -3.0633   13.0180   15.2556   17.4432   19.4943  -16.0663   -2.4733   11.4157   17.4134  -10.6778  -16.9950    5.1923   12.0087
    5.2316    2.0046   14.1176   -5.0475   17.1182   -6.1757    8.4378   11.9649    4.4242   22.8422  -23.0485   -5.9786   23.2909  -20.4804   17.4655    8.6725
    5.2762   -8.4564   18.2362  -10.4512   12.2180   17.1640   22.7231   -7.9217  -16.6471  -11.9341    5.3935   10.1627   -1.7822  -18.4441   10.6187    5.9201
   11.8713   10.3754   16.9336   -8.8719   13.9635    8.5727   12.7539  -14.0447  -28.7104   10.4878    0.0535    2.0597    8.6829   -4.2881   27.7403   16.2610
   -5.2706    3.7313   21.6617  -21.4129   29.6162   11.3101   27.9497    9.8417  -22.4691    8.1512    0.6769   -1.0082    6.4324   -0.3772   26.5080   12.3203
  -25.3620   -5.9851   11.3456  -29.9013   11.8891   18.3878   17.7220    8.6295  -15.0348   11.7951  -34.1147    9.5656  -16.6761  -25.3711   15.9841    2.6499
   32.0527  -14.4586   22.4979    8.2390   37.6135   30.3129   12.6405   12.1462   -1.7126   16.4043   34.6161    9.9245   27.6886   -2.3919   21.9304   22.3448
   18.7039   12.7307    7.6489  -19.3246   19.7053   36.0388   -1.7405   -9.1993  -11.7371   14.9111   42.9560   16.0717    2.6281  -24.3129   12.1685   50.5756
   32.5006   -3.2196   15.6657   -9.4225   34.4549   52.2375    9.7007   22.1439  -12.6050   16.1152   24.9674   17.9730    8.6456   13.1522   32.5255   63.2251
   22.4765   -4.8659   29.0310    7.2803   44.9706   25.2533    7.6699   35.2114  -10.1218    8.6039   50.4770   37.4022   12.2315  -16.0756   36.3739   19.8402
   48.2240   16.2211   28.3039   -0.6217   43.5139   50.5793   -0.6533   11.9925  -11.5865    8.3372   46.0308   32.4310   -8.8339    8.8712   -0.4805   48.0733
   11.7047   -3.9472   27.0932    3.4864   39.6667   24.1539   18.0517   25.3794    5.4614   -9.9581   21.6964   17.3386   27.8180   16.4992   21.5974   42.1414
   -9.9286    6.1899   61.6588   -2.2494   29.2529   17.5899   27.2931   35.5218   36.7595   24.3471    2.9909   29.9213   31.3656   28.6066   14.4134    4.4372
   -1.5634   -2.5133   61.0924    7.2582   15.6025    3.2662   12.7704   65.2843    8.2906    1.7293   -6.3853   29.6491   46.7374   19.6381   20.3943   24.4863
   24.6176   12.4974   52.7100    7.2313   33.7640   25.2432   10.6871   33.1016   14.9972   19.4886    8.1686   36.3052   18.0770   -1.4993   11.0357   10.7673
    6.5605   -5.7408   23.4070  -15.6298   25.2018   16.9284   30.7856   20.4674   -7.2720   15.8910   15.0681   20.8196   12.4573    1.4726   27.4948   15.4116
   11.8954   -4.0401   21.0737  -11.1100   18.6080   27.1161   17.8817   -4.3233   -3.5467   18.9405    9.1683   32.3053   -3.8770  -14.1224   17.9241   16.5994
   12.1514   16.5679   14.2502   -7.4974   37.0718   49.7862    9.8393   26.7070    7.8551    8.9037   -7.0270   24.7438   -1.0138  -20.8358   19.4076   37.4557
    7.7048   15.9947   22.4665   -5.6239   18.2651   21.4577   -7.0290   34.1455   -7.8212   11.1206   10.3417   27.2599    8.6751   -9.7987   39.7746   30.1916
   -6.8981   -9.3803   44.4886  -30.1471   40.0425   23.8180   16.7913    6.1967   -9.7867   19.2595    3.7051   22.4449   -4.2351  -10.8957   17.4246   19.7914
    0.4337    9.1229   23.4034   -9.5021   12.4409   11.0318    3.1281   22.7158  -36.5507   29.7193    4.3535   27.0945    3.5542  -16.2200   37.8575   36.8338
   27.0904   -5.5272   25.3071    4.5881   15.9523   14.2232   11.5042   25.3785    5.4288   30.2172    4.8587   46.3357    2.0642  -36.7948   21.7841   24.2807
   11.9389    3.0703   17.4522  -25.7850   35.8482   22.3387   10.9844   18.5763   -8.4529  -25.6632   14.8003   27.3198    6.5248   -9.0453   25.2100   33.8112
   25.0843    5.3389   42.5562    4.3787   16.0434   33.1410    8.1781   19.7293  -14.8061   15.7858   24.5224   16.0153   25.3372    6.0971    6.7747   56.8271
   15.7529    3.5825   33.8347   -6.9395   37.0054   37.8564   17.3083   30.6840    5.7490   10.7653   31.4631   31.6675    2.7187    8.4376   25.2807   40.9887
   40.0427   16.6011   32.5557    3.6769   28.5854   43.1517   12.8803   27.4759   11.9094   14.9748   34.3143   54.1567    7.9174  -12.3617   32.4838   68.2135
   23.2907   -0.1903   23.1127   -8.1812   33.4546   54.0454   41.5013   22.2376   -0.5338   32.6714   54.7652   18.9165    8.6791   17.8434   41.9372   57.5749
   29.0013    7.3588   62.8455    9.7099   39.4811   55.0241   14.7783   46.0747   29.7710   30.6601   25.3837   32.9911   28.7622   12.8625   32.8878   52.3222
    0.5320   17.5797   57.1094   13.4973   36.5768   24.0481   21.2357   62.5468   21.4502   16.3942    0.8265   27.2666   43.7973   39.6455   34.0867   42.4546
   13.1713    8.8928   80.4872   41.3605   49.2912   12.3685   17.8028   46.1485    7.6889    5.6170   -7.2710   38.9962   40.2377   26.6740   10.1128   35.2287
   22.4834  -20.1894   56.5042   -1.2319   37.0391   40.6007    8.4682   46.8078    9.3974   22.0760   13.9365   28.4277   37.6203   22.8733   33.4697   42.3178
   23.8335  -12.2873   42.3525  -12.8659   22.0899   42.7007   30.3300   25.4185  -26.7167   24.7509   15.7074   42.8802   -4.5030  -18.9023   21.2668   12.8722
   14.3064  -22.3380   24.1042    0.4071   30.3386   33.6773   17.9427   30.6428  -19.2950   21.0647   19.4057   22.8239   30.7307    7.1557   44.4171   42.0231
   20.5359   11.1455   39.4244  -13.8891   27.6906   34.6136   22.2834   31.1930   -7.6862   27.4395    6.9887   -1.3060   15.7327    3.0280   27.1638   27.7310
    0.1172  -10.7539   48.1365   11.0832   33.3089   25.0949    8.6396   23.0315  -21.9054    4.7065   16.5556   21.8430   20.6695    2.0326   30.2804   42.1249
    1.3981  -11.6075   36.6895  -12.3609   20.7054   36.8568   22.5310   51.0005  -14.0872   27.1359   39.9501   15.6667   12.7554  -11.3591   27.1039   30.3997
   23.2718   -1.1503   33.4863    2.8999   32.8178   32.6537    4.9498   41.1391   11.4604   32.1784   13.4726   36.7001   14.0325   -0.0180   62.6150   26.9415
   27.9570    2.6920   43.4339  -18.1471   45.4742   25.2782   29.3600   21.3776   -6.6007   15.9537    6.8576   39.6990   25.3590  -23.6683    4.9698   41.6199
   16.1674   10.0107    4.0620   -3.1911   30.1249   46.8296   21.4719   37.5651    3.1702   11.5687   17.7723   31.8344  -18.9432   12.0868   32.5121   32.7189
   40.8265   10.0837   69.1454   -5.8108   21.7097   56.7448   11.5184   39.6289  -18.4091    4.2910   28.6306   52.5340    3.3894  -10.0632    2.4675   47.1224
   33.8899   -6.9326   27.4529   -1.9221   24.3908   65.6491   17.0124   30.3211  -22.1544   10.5172   45.6400   34.6834   -1.6323   -8.8601   44.5503   91.1406
   63.9031   12.0558   34.7692   24.4283   49.6673   67.9541   16.3649   30.6071  -11.1165   31.2922   53.9543   42.1928   -0.9608   -4.1001   37.6300   82.9257
   35.8485    0.2185   36.2257    0.6254   55.6796   52.7305   10.3208   25.7968    0.0800    4.6125   26.2251   29.3360   32.3145   -3.4877   38.9166   63.7548
   36.5529   12.8773   48.0199    0.9232   66.2201   78.7919   21.2794   59.1893   19.1742   11.7765   45.5641   38.1554   40.4082   26.4788   20.4154   60.9404
   25.8568    8.0317   76.6022   36.8350   37.0220   38.9024   37.4910   53.5126   22.4176    3.8123   21.7685   53.4600   47.6454   20.3109   38.3688   39.0969
   30.7001  -25.7543   77.3110    8.1337   39.9296   52.4349    3.2349   50.6712    8.3374   36.1996   31.8489   39.1083   30.7998   29.3989   36.0193   54.5720
   26.2259   -1.1482   77.4142    7.2508   33.6739   43.5341   36.7071   50.0768    8.7048   28.9675   22.6506   66.7360   53.2672   41.9382   11.6509   54.5165
   22.1079    0.8941   52.0851    8.8913   38.8679   28.5106    3.0630   32.5234    5.3372   16.5439   36.0621   37.9894   28.9366   -4.1390   27.9779   34.9513
   22.1728   14.2901   67.3178   14.9469   16.7645   60.8317   39.6296   45.6104    2.6352    6.8425   22.2077   39.9696   18.7849   15.1160   12.7300   76.8982
   38.6032    6.8664   60.4540   16.1430   53.8365   10.2164    8.4402   37.7016   10.8189   19.4802   -0.0967   51.2623   13.6615   -3.7959   29.7965   48.7979
   62.4589    8.4071   48.0489   -7.6256   44.8543   42.5621   27.7182   27.9053   16.4969   30.3643   48.1805   64.4205   -2.8168   -3.0982   42.4186   48.9524
   26.7562    0.7352   58.2004   -4.7299   26.0813   30.6875   26.6642   38.2456  -19.6434   35.2730   22.8807   31.3073   15.4652   18.4261   56.6563   57.1089
   38.3312   -8.4534   62.9273    9.3807   52.2420   80.8489   28.5975   47.7899  -30.8045   16.1683   34.2330   43.1708  -13.9535  -16.6378   47.5698   77.8496
   14.7782   11.2284   85.2904    3.9041   43.6736   33.1240   13.1411   39.6193  -12.1961   21.0457   26.6608   39.0134    8.5549    2.5965   34.5834   70.4500
   23.0535   -3.6337   51.0256   14.4908   63.5785   58.1665   38.3042   51.2856  -23.5068   47.1376   16.8601   54.8893    0.8416    5.9338   63.2048   64.4326
   47.8066  -20.6788   35.7718    3.8668   47.7542   57.6831    9.0494   53.8643  -20.4661   17.2257   36.2124   52.4701   12.5466  -16.5396   51.1744   75.7739
   61.9288   -1.4475   80.2627   21.0490   22.1388   89.0515   31.9219   57.2534  -12.0193   20.0181   57.9120   72.1285   33.9760    6.6444   48.2800   76.7844
   55.0679   12.9407   86.7077   -2.3175   33.1572   97.8413   43.0413   25.6644   -7.1630   11.9083   52.7288   29.6658    7.0699   27.2360   38.2763  113.1268
   41.9796   25.9244   67.3331   10.9110   34.4759   87.9480   12.4348   64.1040  -17.8036   19.0024   93.7237   41.6857   34.2792    8.3833   53.2754   60.2803
   46.6295   20.8478   72.5256   20.5624   67.3870   73.0537   17.0154   63.6638    2.7862   26.3854   26.9608   52.5438   42.0148    1.4155   49.3604   74.7587
   25.4846   13.1104   72.3767   43.2994   47.6922   59.5966   21.9104   74.1970   40.3135   41.6451   12.4384   54.8609   60.4085   33.2939   45.0708   67.0835
   44.1731   16.8450   72.9217   20.1836   62.2134   51.4011   31.6578   86.2918   16.8398    1.8710   39.2492   57.5370   58.7117   14.4091   65.0713   75.6083
   37.0334   -1.0370   81.7761   -6.6982   65.1467   56.7019   30.0256   59.6978   26.5976   25.1796   52.6572   52.1576   46.3367   25.4472   67.6853   66.6005
   64.1511   29.4505   63.9094   -1.5457   49.1426   44.5281   23.4791   55.7447  -11.5464   30.1203   41.3939   43.5556   16.6802   -6.0822   43.2975   82.7599
   21.6955    2.6728   86.9173  -22.6883   49.0831   70.6480   30.8610   32.1743  -33.0313    3.2668   20.7686   34.7866    6.8248   -9.1948   50.3439   60.7059
   20.2942    0.9122   81.9206    0.4120   52.2536   60.2528   40.5635   47.2591   30.6423   25.6978   48.4873   67.1468   12.6320    8.7226   28.4707   77.3665
   15.8225  -16.1372   49.4029   11.0338   73.8911   70.3225   34.1949   49.3245    6.6248    0.0989   35.6594   58.0252   17.7109  -31.6470   38.0403   77.6822
   13.4326    5.5421   75.9677   -0.9630   36.1790   54.7340   33.8951   10.5596    7.9834   29.9910   38.1025   26.5381   29.9075   14.5821   24.7835   34.8847
   55.2570   -0.6104   67.4944    1.0358   58.2666   74.5739   14.4102   54.8437   -1.4807    9.5720   34.5253   72.2491   10.9814    3.7234   28.6942   64.2097
   25.8600  -11.2847   65.8732   -2.2418   30.4463  102.5025   38.9948   47.8066    5.6584   31.2944   58.9619   50.1075   -0.7704   -6.5752   31.7702   47.3991
   41.6911  -37.9499   81.7125   16.7606   65.7901   70.2983   26.5532   62.6128   15.4922   46.0143   51.3473   50.7863   21.3817   -2.3776   60.6523   99.1720
   55.4562   -4.1135   40.9375   -1.1727   86.7854  110.9253   44.2076   46.4384   10.6109   22.4377   51.1509   32.6764   30.2924   26.0994   62.9320   91.0542
   49.5487   -0.0754   79.0434   16.0336   57.1886  112.5752   31.6842   75.2740    2.4734    1.9589   77.0983   48.6567   28.8413  -16.4876   47.6441   96.2642
   75.6847    0.9430   62.3475   21.2493   47.9612   94.0197   26.6917   34.9824  -11.8800   32.2776   85.9274   75.6704   49.2666    4.2794   63.6199   97.5433
   49.4726    7.1422   72.4722   16.2629   58.3611   93.3372   42.5572   61.5126   26.8121   30.2621   60.8382   29.2638   15.9643   -1.3031   23.5744  111.0243
   52.3909   -0.1736  118.1940   30.8980   58.8924   92.2075   48.3393   69.2698    2.8613   26.1032   59.1911   54.2439   36.5362  -11.7682   37.0450   67.0651
   29.7631  -26.4196  109.3549   30.8910   67.5949   63.4849   44.3775   79.7190   42.0596   52.9128   47.7569   92.1125   57.5256   32.4243   66.1160   68.0977
   27.8746   28.3217   98.0203   50.1700   61.0438   62.9252   30.8254   78.9920   14.2901   32.6159   51.7623   63.6022   58.2691   13.2113   60.9376   83.6374
   30.4841   14.5326   98.7305   25.5066   49.9342   85.8722   40.9879   70.3123   14.2762   30.7618   46.9186   70.6985   43.5227   14.2773   45.1636   85.8412
   52.9206   23.4607   72.7510   -2.3903   66.7806   73.8898   18.6056   76.6227  -12.0419   27.1603   47.2126   80.4097   -0.8934   17.2591   40.7144   65.6636
   42.4912   -3.7901   81.8765    6.1406   51.0845   79.4289   37.8950   54.1332  -16.7364   28.7075   52.0742   51.2165   27.9005   -5.7362   37.0402   79.7461
   68.8148  -11.1616   85.0842   18.0636   52.6744   80.6601   47.8998   56.8996   11.3549   45.3050   51.7563   56.3160   15.2280   33.9184   48.8860   76.3492
   26.1859    6.4818   64.3392   10.1072   61.9956   66.8068   32.7720   53.3708  -25.5721   36.1160   48.6967   59.5169   16.7174    0.7266   65.1207   61.1794
   41.5048   -7.1848  107.7531   15.2535   63.1957   79.7791   29.1948   81.1423  -24.0731   50.8461   66.4005   80.2223   41.0844    2.2589   56.2460   72.7431
   62.1953   -2.3312   85.1567    0.7820   62.7112   53.1003   31.8949   61.9279    9.9182   19.3671   27.1182   80.0963   58.2098   11.4245   74.6251   88.0028
   54.0356   25.7917   85.2338   23.7979   79.3182   61.2895   37.4489   75.2051  -12.5958   50.2469   56.8445   60.3289   30.9387    2.0235   83.5576  121.5881
   53.8628    6.4370   81.5390   -9.6815   82.3268   75.5936   20.9419   62.9699   -0.4838   40.7315   53.3830   81.3871   36.8807   19.8819   59.9470  141.3720
   67.9088    6.8290   77.0308   25.5577   76.8479  106.4009   14.4749   62.6891    1.8040   13.5469   41.7770   51.3083   31.1931   19.1352   60.9333  103.1614
   49.3969    2.8244   93.3059   17.9222   73.1217  126.1972   39.8899   81.5314  -23.0058   36.6587   81.8170   76.0654   42.9766   24.3119   77.8808  123.6808
   87.5701  -10.6321   81.8787   16.3714   61.4169  110.0872   17.5287   65.6212  -14.3070   17.3709   92.9510   82.6920   52.4376   -1.8078   70.7357   77.7866
   69.1134    9.1870   74.0543  -32.9913   48.3059  133.1157   15.4133   69.2022  -13.4041   31.7301   51.2093   76.2619   40.6770    1.0643   70.3604  126.3186
   74.6844   14.3529   87.3076   36.1650   56.1860  106.0070   53.0831  105.8045  -12.7302   37.3699   35.4565   81.7094   36.6510   18.7082   40.8359  109.3452
   43.7092   18.4087  100.9821   27.8748   70.2665   90.1227   51.7600   79.2198   45.1923   74.8707   47.1076   63.5877   79.3959   59.0766   71.1994   58.6034
   37.9578   12.7600  115.0494    5.8827   92.4287   68.2695   24.4935   87.3104   40.4750   19.6072   24.2387   72.6732   45.5105   50.7274   61.4893  101.5336
   58.3072   19.3893  119.8139   11.3014   74.9190  108.0919   22.5453   75.5394   19.4822   37.9617   59.1035   86.7762   37.9443   23.8556   60.1143   89.0370
   45.2518    8.3525   91.6176   51.2865   73.8824  101.1397   20.7842   51.5311   -2.3402   40.7902   54.6768   62.7699   16.8668    1.5064   68.6067  101.5635
   31.2689  -26.3611   62.6503    4.0199   42.5012   64.4562   22.1972   56.1730   10.9980   31.5066   22.8467   55.8045   59.5530    1.0387   60.3656   98.8959
   46.8130    9.1659   54.1549   -3.6252   79.3574   63.2073   27.4336   74.1022  -28.0544   31.4303   55.8800   38.9311   32.4726   39.6569   67.7920  119.9636
   51.3537    5.1194  117.2743   17.8514   44.7308   91.8493   19.6796   77.6647   17.0059   31.9182   59.5845   64.9742   18.6273  -16.0378   52.3926   93.9586
   48.1280   17.6353   73.4424   25.8187   94.4469   67.9910   40.7395   90.1688   24.7098   18.1938   50.7534   71.1009   72.2773    0.5075   82.0469  104.1610
   47.2697  -17.6449   94.2619   -3.0961   59.3448   73.6039   39.0206   63.4725  -16.7944   -0.3748   36.8124   64.4457   36.8232    3.8152   63.3479   77.4674
   17.2559   14.7662  110.3816   25.8861   79.4889   68.1699   20.8046   64.3010   -4.8506   29.2133   61.3732   33.3051   19.7864   14.6357   80.4661   84.8510

% Copyright (c) 2000, Daniel B. Rowe
% Performs the vectorization operation
% Takes a matrix and turns it into a vector
% by stacking its columns.

function [ vector ] = vec(matrix) 
[p,m]=size(matrix);
vector=zeros(p*m,1);
for cols=1:m
for rows=1:p
vector(rows+p*(cols-1),1 )=matrix(rows,cols); 
 end
end

% Copyright (c) 2000, Daniel B. Rowe
% Performs the trace operation
% Adds the diagonal elements of a square matrix.

function [ trace ] = tr(matrix) 
[p,p]=size(matrix);
trace=0;
for diagelemt=1:p
 trace=trace+matrix(diagelemt,diagelemt); 
end

-

-/2

/2