function [LogMargs,CovProb] = MargLikeModifiedHarmonic(PostSample,LogPost,ComputeSequential,CovProb,DSGEModel,WaitHandle)
% MargLikeModifiedHarmonic: Computes the marginal likelihood from a MCMC sample using
%                           the Modified Harmonic Mean estimator.
%
% USAGE:
%
%       [LogMargs,CovProb] = MargLikeModifiedHarmonic(PostSample,LogPost,ComputeSequential,CovProb,DSGEModel,WaitHandle)
%
% REQUIRED INPUT:  PostSample (matrix) of dimension (NumIter x NumParam) with posterior draws of parameters.
%
%                  LogPost (vector) of dimension (NumIter x 1) with log posterior values for the parameters.
%
%                  ComputeSequential (boolean), which determines if the marginal likelihood should be calculated
%                              sequentially or just for the full sample
%
%                  CovProb (vector), with coverage probabilities.
%
%                  DSGEModel (structure), whose fields contain all required information
%                              about the model, such as input files, variable and sample
%                              data. Actual observations are not stored in this file.
%
% OPTIONAL INPUT:  WaitHandle (handle) to a wait dialog. Allows the user to end the function prematurely.
%
% REQUIRED OUTPUT: LogMargs (matrix) with all calculated marginal likelihoods. The values
%                              for the coverage probabilities are located in the columns
%                              and for the different samples in the rows. The sample size
%                              for the sequential estimate is located in the first column.
%
%                  CovProb (vector), with coverage probabilities.
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2008 European Central Bank.
%
%                       First version: August 30, 2006.
%                        This version: December 19, 2008.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% LICENSE INFORMATION:
%
%      YADA is free software: you can redistribute it and/or modify
%      it under the terms of the GNU General Public License as published by
%      the Free Software Foundation, either version 3 of the License, or
%      (at your option) any later version.
%
%      This program is distributed in the hope that it will be useful,
%      but WITHOUT ANY WARRANTY; without even the implied warranty of
%      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%      GNU General Public License for more details.
%
%      You should have received a copy of the GNU General Public License
%      along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%      YADA is released under the GNU General Public License, Version 3,
%      29 June 2007 <http://www.gnu.org/licenses/>. The current release of
%      the program was last modified by the ECB on the "This version" date
%      above.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANGELOG:
%
% * 31-8-2006: Added everything below line 65.
%
% * 1-9-2006:  Updated the documentation.
%
% * 9-11-2006: Took the name change chi2inv to YADAchi2inv into account.
%
% * 29-11-2006: Fixed some output bugs with LogMargs.
%
% * 12-12-2006: Added a drawnow command at the end of the sequential estimation loop.
%
% * 14-12-2006: Added the cancel call check.
%
% * 21-12-2006: Added a missing N to the expression "LogMargs = [LogMargs;N NaN*ones(1,length(CovProb))];".
%
% * 8-1-2007:   If "SeqGrid" is empty, it is now set to NumIter.
%
% * 8-2-2007:   Fixed a bug related to clicking on the cancel button. Missing "StopEstimation" test for storing
%               data.
%
% * 13-11-2007: Updated the documentation.
%
% * 27-3-2008:  Made sure that the sequential grid has at least 2 values.
%
% * 23-5-2008:  Updated the documentation.
%
% * 21-7-2008:  Added the "UseYADAVersion" internal variable.
%
% * 19-12-2008: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

if (exist('WaitHandle')==0)|(nargin==5);
   WaitHandle = [];
else;
   WaitControls = get(get(WaitHandle,'Parent'),'UserData');
   WaitTxt = get(WaitControls.text,'String');
   k = size(WaitTxt,1);
   kLine = WaitTxt{k};
end;
%
if isempty(CovProb)==1;
   %
   % Compute coverage probabilities based on the DSGEModel structure
   %
   CovProb = (DSGEModel.CovStartValue:DSGEModel.CovIncValue:DSGEModel.CovEndValue);
end;
CovProb = vec(CovProb)';
%
[NumIter,NumParm] = size(PostSample);
%
if ComputeSequential==1;
   Origin = DSGEModel.SequentialStartIterationValue;
   %
   % make sure that Origin is less than NumIter
   %
   if Origin==NumIter;
      Origin = floor(NumIter/10);
   end;
   StepLength = DSGEModel.SequentialStepLengthValue;
else;
   Origin = NumIter;
   StepLength = 1;
end;
%
% Create the sequential grid over which we will compute the marginal likelihood
%
SeqGrid = (Origin:StepLength:NumIter);
if isempty(SeqGrid)==1;
   SeqGrid = NumIter;
elseif SeqGrid(length(SeqGrid))<NumIter;
   SeqGrid = [SeqGrid NumIter];
end;
%
% check how critical values from the chi-square should be computed
%
if exist('YADAchi2inv')==2;
   UseYADAVersion = 2;
else;
   if exist('chi2inv')==2;
      UseYADAVersion = 1;
   else;
      UseYADAVersion = 0;
   end;
end;
%
% initialization
%
SumSq = zeros(NumParm,NumParm);
SumP = zeros(NumParm,1);
Sigma = zeros(NumParm,NumParm);
MaxLogPost = max(LogPost);
mu = zeros(NumParm,1);
StartEst = 1;
LogMargs = [];
StopEstimation = 0;
for N=SeqGrid;
   %
   % update the wait dialog
   %
   if isempty(WaitHandle)==0;
      WaitTxt{k} = [kLine ' Sample size: ' int2str(N) ' (' int2str(NumIter) ')'];
      set(WaitControls.text,'String',WaitTxt);
   end;
   %
   % calculate the mean and the coveriance matrix
   %
   ValidSigma = 1;
   SumP = SumP+sum(PostSample(StartEst:N,:))';
   mu = SumP/N;
   SumSq = SumSq+PostSample(StartEst:N,:)'*PostSample(StartEst:N,:);
   Sigma = (SumSq/N)-(mu*mu');
   %
   % Check if Sigma is valid
   %
   SigmaTest = sum(log(eig(Sigma)));
   if (isinf(SigmaTest)==1)|(isnan(SigmaTest)==1)|(isreal(SigmaTest)==0);
      LogMargs = [LogMargs;N NaN*ones(1,length(CovProb))];
      ValidSigma = 0;
   end;
   %
   % We can now try to calculate the marginal likelihood
   %
   if ValidSigma==1;
      InvSigma = inv(Sigma);
      LogMarg = [];
      %
      % we look at one coverage probability at a time
      %
      for p=CovProb;
         if UseYADAVersion==2;
            CritValue = YADAchi2inv(p,NumParm);
         elseif UseYADAVersion==1;
            CritValue = chi2inv(p,NumParm);
         else;
            CritValue = YADACheapGammaInv(p,(NumParm/2))*2;
         end;
         TempValue = 0;
         for i=1:N;
            if isempty(WaitHandle)==0;
               %
               % check if the cancel button was pressed
               %
               if strcmp(get(WaitHandle,'UserData'),'cancel')==1;
                  StopEstimation = 1;
                  break;
               end;
            end;
            Deviation = (PostSample(i,:)-mu')*InvSigma*(PostSample(i,:)'-mu);
            if Deviation<=CritValue;
               LogTruncNormalDens = -log(p)-(((NumParm*log(2*pi))+SigmaTest+Deviation)/2);
               TempValue = TempValue+exp(LogTruncNormalDens-LogPost(i)+MaxLogPost);
            end;
         end;
         LogMarg = [LogMarg (MaxLogPost-log(TempValue/N))];
         drawnow;
         if StopEstimation==1;
            break;
         end;
      end;
      if StopEstimation==0;
         LogMargs = [LogMargs;N LogMarg];
      end;
   end;
   StartEst = N+1;
   if StopEstimation==1;
      LogMargs = [];
      break;
   end;
end;

%
% end of MargLikeModifiedHarmonic.m
%
