wolffd@0: function [B,B2,dist] = parzen(data, mu, Sigma, N)
wolffd@0: % EVAL_PDF_COND_PARZEN Evaluate the pdf of a conditional Parzen window
wolffd@0: % function B = eval_pdf_cond_parzen(data, mu, Sigma, N)
wolffd@0: %
wolffd@0: % B(q,t) = Pr(data(:,t) | Q=q) = sum_{m=1}^{N(q)} w(m,q)*K(data(:,t) - mu(:,m,q); sigma)
wolffd@0: % where K() is a Gaussian kernel with spherical variance sigma,
wolffd@0: % and w(m,q) = 1/N(q) if m<=N(q) and = 0 otherwise
wolffd@0: % where N(q) is the number of mxiture components for q 
wolffd@0: %
wolffd@0: % B2(m,q,t) =  K(data(:,t) - mu(:,m,q); sigma) for m=1:max(N)
wolffd@0: 
wolffd@0: % This is like eval_pdf_cond_parzen, except mu is mu(:,m,q) instead of mu(:,q,m)
wolffd@0: % and we use 1/N(q) instead of mixmat(q,m)
wolffd@0: 
wolffd@0: if nargout >= 2
wolffd@0:   keep_B2 = 1;
wolffd@0: else
wolffd@0:   keep_B2 = 0;
wolffd@0: end
wolffd@0: 
wolffd@0: if nargout >= 3
wolffd@0:   keep_dist = 1;
wolffd@0: else
wolffd@0:   keep_dist = 0;
wolffd@0: end
wolffd@0: 
wolffd@0: [d M Q] = size(mu);
wolffd@0: [d T] = size(data);
wolffd@0: 
wolffd@0: M = max(N(:));
wolffd@0: 
wolffd@0: B = zeros(Q,T);
wolffd@0: const1 = (2*pi*Sigma)^(-d/2);
wolffd@0: const2 = -(1/(2*Sigma));
wolffd@0: if T*Q*M>20000000 % not enough memory to call sqdist
wolffd@0:   disp('eval parzen for loop')
wolffd@0:   if keep_dist,
wolffd@0:     dist = zeros(M,Q,T);
wolffd@0:   end
wolffd@0:   if keep_B2
wolffd@0:     B2 = zeros(M,Q,T);
wolffd@0:   end
wolffd@0:   for q=1:Q
wolffd@0:     D = sqdist(mu(:,1:N(q),q), data); % D(m,t)
wolffd@0:     if keep_dist
wolffd@0:       dist(:,q,:) = D;
wolffd@0:     end
wolffd@0:     tmp = const1 * exp(const2*D);
wolffd@0:     if keep_B2,
wolffd@0:       B2(:,q,:) = tmp;
wolffd@0:     end
wolffd@0:     if N(q) > 0
wolffd@0:       %B(q,:) = (1/N(q)) * const1 * sum(exp(const2*D), 2);
wolffd@0:       B(q,:) = (1/N(q)) * sum(tmp,1);
wolffd@0:     end
wolffd@0:   end
wolffd@0: else
wolffd@0:   %disp('eval parzen vectorized')
wolffd@0:   dist = sqdist(reshape(mu(:,1:M,:), [d M*Q]), data); % D(mq,t)
wolffd@0:   dist = reshape(dist, [M Q T]);
wolffd@0:   B2 = const1 * exp(const2*dist); % B2(m,q,t)
wolffd@0:   if ~keep_dist
wolffd@0:     clear dist
wolffd@0:   end
wolffd@0:   
wolffd@0:   % weights(m,q) is the weight of mixture component m for q  
wolffd@0:   %    = 1/N(q) if m<=N(q) and = 0 otherwise
wolffd@0:   % e.g., N = [2   3   1], M = 3,
wolffd@0:   % weights = [1/2 1/3 1   = 1/2 1/3 1/1      2 3 1     1 1 1
wolffd@0:   %            1/2 1/3 0     1/2 1/3 1/1 .*   2 3 1 <=  2 2 2
wolffd@0:   %            0   1/3 0]    1/2 1/3 1/1      2 3 1     3 3 3
wolffd@0:    
wolffd@0:   Ns = repmat(N(:)', [M 1]);
wolffd@0:   ramp = 1:M;
wolffd@0:   ramp = repmat(ramp(:), [1 Q]);
wolffd@0:   n = N + (N==0); % avoid 1/0 by replacing with 0* 1/1m where 0 comes from mask
wolffd@0:   N1 = repmat(1 ./ n(:)', [M 1]);
wolffd@0:   mask = (ramp <= Ns);
wolffd@0:   weights = N1 .* mask;
wolffd@0:   B2 = B2 .* repmat(mask, [1 1 T]);
wolffd@0:   
wolffd@0:   % B(q,t) = sum_m B2(m,q,t) * P(m|q) = sum_m B2(m,q,t) * weights(m,q)
wolffd@0:   B = squeeze(sum(B2 .* repmat(weights, [1 1 T]), 1)); 
wolffd@0:   B = reshape(B, [Q T]); % undo effect of squeeze in case Q = 1
wolffd@0: end
wolffd@0: 
wolffd@0:   
wolffd@0: