% Class for normalian Wishart distribution used by Dirichlet Process with
% SV.

% The novelty is we keep track of U = chol(S+kappa*(mu*mu')) instead of S.
% S is the H matrix in equation (B.73) in Zellner (1971, p.396) which is
% the precision matrix in standard form.

% Besides the Data inputs (X), we also need the SV component inputs (h).

% Univaraite case d = 1.
classdef normalWishartSV
        properties
         kappa_
         mu_
         nu_
         U_ % we track U, instead of S, the precision matrix
     end
     
     methods
         function obj = normalWishartSV(kappa,mu,nu,S) % hyperparameters
             U = chol(S+kappa*(mu*mu'));
             obj.kappa_ = kappa;
             obj.mu_ = mu;
             obj.nu_ = nu;
             obj.U_ = U;
         end
         
         function obj = clone(obj)
         end
         
         function obj = addData(obj, X, h) % X is d by n, h is 1 by n
             kappa0 = obj.kappa_;
             mu0 = obj.mu_;
             nu0 = obj.nu_;
             U0 = obj.U_;
             
             exph = exp(-h/2);
             exph2 = exp(-h);
             n = size(X,2);
             
             kappa = kappa0+sum(exph2);
             mu = (kappa0*mu0+sum(exph2.*X,2))/kappa;
             nu = nu0+n;
             U = chol(U0'*U0+(exph.*X)*(exph.*X)');
             
             obj.kappa_ = kappa;
             obj.mu_ = mu;
             obj.nu_ = nu;
             obj.U_ = U;
         end
        
         function obj = addSample(obj, x, h) %x is d by 1
             kappa = obj.kappa_;
             mu = obj.mu_;
             nu = obj.nu_;
             U = obj.U_;
             
             exph = exp(-h/2);
             exph2 = exp(-h);
             
             kappa = kappa+exph2;
             mu = mu+exph2*(x-mu)/kappa;
             nu = nu+1;
             U = cholupdate(U,exph*x,'+');
             
             obj.kappa_ = kappa;
             obj.mu_ = mu;
             obj.nu_ = nu;
             obj.U_ = U;
         end
         
         function obj = delSample(obj, x, h)
             kappa = obj.kappa_;
             mu = obj.mu_;
             nu = obj.nu_;
             U = obj.U_;
             
             exph = exp(-h/2);
             exph2 = exp(-h);
             
             kappa = kappa-exph2;
             mu = mu-exph2*(x-mu)/kappa;
             nu = nu-1;
             U = cholupdate(U,exph*x,'-');
             
             obj.kappa_ = kappa;
             obj.mu_ = mu;
             obj.nu_ = nu;
             obj.U_ = U;
         end
         
         function y = logPredPdf(obj, X, h) 
             % predictive densiity is multivariate T(mu, S*(Kappa+exph2)/(Kappa*v)
             % v is degrees of freedom = nun-d+1 
             kappa = obj.kappa_;
             mu = obj.mu_;
             nu = obj.nu_;
             U = obj.U_;
             
             exph = exp(-h/2);
             exph2 = exp(-h);
             d = size(X,1);
             v = (nu-d+1);
             U = sqrt((1+exph2/kappa)/v)*cholupdate(U,sqrt(kappa)*mu,'-');
             
             X = bsxfun(@minus,exph.*X,exph.*mu);
             Q = U'\X;
             q = dot(Q,Q,1);  % quadratic term (M distance)
             o = -log(1+q/v)*((v+d)/2);
             c = gammaln((v+d)/2)-gammaln(v/2)-(d*log(v*pi)+2*sum(log(diag(U))))/2;
             y = c+o;
         end
     end
end
