%%%% Code to compute the kernel between two Gaussian emission HMMs function [ba] = bhatt(mu1,mu2,cov1,cov2) % function [ba] = bhatt(mu1,mu2,cov1,cov2) % Calculates the bhattacharya affinity for two distributions % mu1, mu2, cov1 and cov2 are the means and covariances % mu1 and mu2 are column vectors, cov1 and cov2 are square % matrices % % Yingbo Song (yingbo@cs.columbia.edu) [n m] = size(cov1); % add a small regularizer to the diagonal, ensure covm is full rank cov1 = cov1 + 1e-5*trace(cov1)*eye(n); cov2 = cov2 + 1e-5*trace(cov2)*eye(n); iC1 = inv(cov1); iC2 = inv(cov2); Cd = inv(iC1 + iC2); Md = iC1*mu1 + iC2*mu2; % calculate the affinity rho = 0.5; % exponential, should be 0.5 for Bhattacharyya affinity normalizer = (2*pi)^((1-2*rho)*(n/2)) * rho^(-n/2); normalizer = normalizer * det(cov1)^(-rho/2) * det(cov2)^(-rho/2) * sqrt(det(Cd)); ba = normalizer * exp((-rho/2)*(mu1'*(iC1)*mu1 + mu2'*(iC2)*mu2 - Md'*(Cd)*Md)); % end of function function [K] = elkernel(hmm1,hmm2,T) % function [K] = elkernel(hmm1,hmm2,T) % Tony's iterative PPK % NOT in logspace if ( sum(hmm1.prior) < 0 ) error('elkernel: The hmm parameters cannot be in log-space!'); end prior1 = hmm1.prior; prior2 = hmm2.prior; transmat1 = hmm1.trans; transmat2 = hmm2.trans; M = length(prior1); N = length(prior2); pad = 0.45; pot = zeros(M,N); for i=1:M for j=1:N pot(i,j) = bhatt(hmm1.model(i).mu,hmm2.model(j).mu,... hmm1.model(i).covv + pad,hmm2.model(j).covv + pad); end end K = 0; if (T==1) for i=1:M for j=1:N K=K+prior1(i)*prior2(j)*pot(i,j); end end else sep1 = zeros(M,N); for i=1:M for j=1:N sep1 = sep1 + prior1(i)*prior2(j)*pot(i,j)*transmat1(i,:)'*transmat2(j,:); end end for t=2:T sep2 = zeros(M,N); for i=1:M for j=1:N sep2 = sep2 + sep1(i,j)*pot(i,j)*transmat1(i,:)'*transmat2(j,:); end end sep1 = sep2; end K = sum(sum(sep2.*pot)); end