diff toolboxes/FullBNT-1.0.7/KPMstats/parzen.m @ 0:e9a9cd732c1e tip

first hg version after svn

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