% PoissonSeriesRep
function [muX sigX] = PoissonSeriesRep(alpha, c, mu, muW, sigW, N)
% Following Lemke, Riabiz and Godsill  (Digital Signal Processing 2015)
% requires alpha > 1.  alpha = 2 may not work.
% N = # of kernels
% muX, sigX ~ 1 x N give means and sds of kernels
% Use at least N = 1000 for good fit.
% Code illustrated and graphed in file PoissonSeriesRepDemo.m
 
% Pre-compute muW, sigW from alpha, beta with code below:
% Calpha = (1-alpha)/(gamma(2-alpha)*cos(pi*alpha/2));
% if abs(beta) == 1
%   muW = sign(beta) * Calpha^(1/alpha)
%   sigW = 0
%   % beta = +/-1 gives OK cdf, but pdf undefined.  Kernels become particles.
% else
%   % Invert FindBetaFromMu to find muStar from beta:
%   betafun = @(muStar,alpha,beta) FindBetaFromMu(muStar,alpha) - beta;
%   bfun = @(muStar) betafun(muStar, alpha, beta);
%   muStar = fzero(bfun,0)
%   EWStarAlpha = GaussRomberg(@(x) abs(x).^alpha, muStar, 1, 1000)
%   sigW = (Calpha/EWStarAlpha)^(1/alpha)
%   muW = sigW*muStar
% end
g = 100; % Poisson truncation point per LRG
muR = [alpha/(1-alpha)*g^((alpha-1)/alpha); ... 
       alpha/(2-alpha)*g^((alpha-2)/alpha)];
SigR = [muR(2)  alpha/(3-alpha)*g^((alpha-3)/alpha); ... 
          0       alpha/(4-alpha)*g^((alpha-4)/alpha)];
SigR(2,1) = SigR(1,2);

% Vectorized version:
Gam = cumsum(random('exp', 1, N, 2*g), 2);    % ~ N x 2g
GamGTg = find(Gam > g);
GamIa = Gam.^(-1/alpha);
GamIa(GamGTg) = 0;
m = sum(GamIa,2);
s = sum(GamIa.^2,2); 
R = mvnrnd(repmat(muR', N, 1), repmat(SigR, 1, 1, N));   % ~ N x 2
muX = c * muW * (m + R(:,1))' + mu;
sigX = c * sigW * sqrt(s + R(:,2))';   % ~ 1 x N

% Loopy version:
% muX =  nan(1,N);
% sigX = nan(1,N);
% R1muW = nan(1,N);
% R2sigW2 = nan(1,N);
% mmuW = nan(1,N);
% ssigW2 = nan(1,N);
% for i = 1:N
%   Gam = random('exp', 1, 2*g, 1);
%   Gam = cumsum(Gam);
%   GamGTg = Gam > g;
%   M = find(GamGTg, 1)-1;
%   Gam = Gam(1:M);
%   m = sum(Gam .^ (-1/alpha));
%   s = sum(Gam .^ (-2/alpha));
%   mmuW(i) = m*muW;
%   ssigW2(i) = s*sigW^2;
%   R = mvnrnd(muR, SigR);
%   R1muW(i) = R(1)*muW;
%   R2sigW2(i) = R(2)*sigW^2;
%   muX(i) = c*(m + R(1))*muW + mu;
%   sigX(i) = c* sqrt(s + R(2))*sigW;
% end
return



  