%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Analytical merge of 2 single-Gaussian processes
% mu: mean
% sigma: std dev
% n: number of samples

function [mu12, sigma12, n12, diag_sigma12] = merge_single_gaussian( mu1, sigma1, n1, mu2, sigma2, n2, unbiased )
  
  if nargin < 6
    error( [ mfilename '.merge_single_gaussian() needs at least 6 input arguments.' ] );
  end
  
  if ~exist( 'unbiased', 'var' )
    unbiased = 1;
  end

  mu12 = ( n1 * mu1 + n2 * mu2 ) / ( n1 + n2 );
  
  n12 = n1 + n2;

  if unbiased
    
    sigma12 = ( ( n1 - 1 ) * sigma1 + n1 * mu1(:) * mu1(:).' + ( n2 - 1 ) * sigma2 + n2 * mu2(:) * mu2(:).' - ( n1 + n2 ) * mu12(:) * mu12(:).' ) / ( n1 + n2 - 1 );
      
  else

    sigma12 = ( n1 * sigma1 + n1 * mu1(:) * mu1(:).' + n2 * sigma2 + n2 * mu2(:) * mu2(:).' - ( n1 + n2 ) * mu12(:) * mu12(:).' ) / ( n1 + n2 );
    
  end
  
  % Make sure that it is positive definite ( for the inversion in log_mvnpdf.m )
  [ R, p ] = chol( sigma12 );

  diag_sigma12 = ( p~= 0 );

  if diag_sigma12
    % If not positive definite, fall back on the diagonal covariance matrix
    sigma12 = diag( diag( sigma12 ) );
  end

  mu12     = full( mu12 );
  sigma12  = full( sigma12 );