diff options
Diffstat (limited to 'sourcecodes/bnt-master/KPMstats')
81 files changed, 4072 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/KPMstats/CVS/Entries b/sourcecodes/bnt-master/KPMstats/CVS/Entries new file mode 100644 index 00000000..9e8985ba --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/CVS/Entries @@ -0,0 +1,80 @@ +/KLgauss.m/1.1.1.1/Tue Apr 26 02:29:16 2005// +/README.txt/1.1.1.1/Tue Apr 26 02:29:18 2005// +/beta_sample.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/chisquared_histo.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/chisquared_prob.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/chisquared_readme.txt/1.1.1.1/Tue Apr 26 02:29:18 2005// +/chisquared_table.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/clg_Mstep.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/clg_Mstep_simple.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/clg_prob.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/condGaussToJoint.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/cond_indep_fisher_z.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/condgaussTrainObserved.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/condgauss_sample.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/convertBinaryLabels.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/cwr_demo.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/cwr_em.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/cwr_predict.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/cwr_prob.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/cwr_readme.txt/1.1.1.1/Tue Apr 26 02:29:18 2005// +/cwr_test.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/dirichlet_sample.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/dirichletpdf.m/1.1.1.1/Sun May 22 23:32:18 2005// +/dirichletrnd.m/1.1.1.1/Sun May 22 23:32:12 2005// +/distchck.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/eigdec.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/est_transmat.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/fit_paritioned_model_testfn.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/fit_partitioned_model.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/gamma_sample.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/gaussian_prob.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/gaussian_sample.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/histCmpChi2.m/1.1.1.1/Tue May 3 20:18:16 2005// +/linear_regression.m/1.1.1.1/Tue Apr 26 02:29:18 2005// +/logist2.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/logist2Apply.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/logist2ApplyRegularized.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/logist2Fit.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/logist2FitRegularized.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/logistK.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/logistK_eval.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/marginalize_gaussian.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/matrix_T_pdf.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/matrix_normal_pdf.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/mc_stat_distrib.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/mixgauss_Mstep.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/mixgauss_classifier_apply.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/mixgauss_classifier_train.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/mixgauss_em.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/mixgauss_init.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/mixgauss_prob.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/mixgauss_prob_test.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/mixgauss_sample.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/mkPolyFvec.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/mk_unit_norm.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/multinomial_prob.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/multinomial_sample.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/multipdf.m/1.1.1.1/Sun May 22 23:32:42 2005// +/multirnd.m/1.1.1.1/Sun May 22 23:32:38 2005// +/normal_coef.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/partial_corr_coef.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/parzen.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/parzenC.c/1.1.1.1/Tue Apr 26 02:29:20 2005// +/parzenC.dll/1.1.1.1/Tue Apr 26 02:29:20 2005// +/parzenC.mexglx/1.1.1.1/Tue Apr 26 02:29:20 2005// +/parzenC_test.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/parzen_fit_select_unif.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/pca.m/1.1.1.1/Tue Apr 26 02:29:20 2005// +/rndcheck.m/1.1.1.1/Tue Apr 26 02:29:22 2005// +/sample.m/1.1.1.1/Tue Apr 26 02:29:22 2005// +/sample_discrete.m/1.1.1.1/Tue Apr 26 02:29:22 2005// +/sample_gaussian.m/1.1.1.1/Tue Apr 26 02:29:22 2005// +/standardize.m/1.1.1.1/Wed May 4 04:35:36 2005// +/student_t_logprob.m/1.1.1.1/Tue Apr 26 02:29:22 2005// +/student_t_prob.m/1.1.1.1/Tue Apr 26 02:29:22 2005// +/test_dir.m/1.1.1.1/Sun May 22 23:32:20 2005// +/unidrndKPM.m/1.1.1.1/Tue May 31 18:19:24 2005// +/unif_discrete_sample.m/1.1.1.1/Tue Apr 26 02:29:22 2005// +/weightedRegression.m/1.1.1.1/Tue Apr 26 02:29:22 2005// +D diff --git a/sourcecodes/bnt-master/KPMstats/CVS/Repository b/sourcecodes/bnt-master/KPMstats/CVS/Repository new file mode 100644 index 00000000..55a108ab --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/CVS/Repository @@ -0,0 +1 @@ +FullBNT/KPMstats diff --git a/sourcecodes/bnt-master/KPMstats/CVS/Root b/sourcecodes/bnt-master/KPMstats/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/KPMstats/KLgauss.m b/sourcecodes/bnt-master/KPMstats/KLgauss.m new file mode 100644 index 00000000..85f00619 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/KLgauss.m @@ -0,0 +1,9 @@ +function kl = KLgauss(P, Q) +%The following computes D(P||Q), the KL divergence between two zero-mean +%Gaussians with covariance P and Q: +% klDiv = -0.5*(log(det(P*inv(Q))) + trace(eye(N)-P*inv(Q))); + +R = P*inv(Q); +kl = -0.5*(log(det(R))) + trace(eye(length(P))-R); + +%To get MI, just set P=cov(X,Y) and Q=blockdiag(cov(X),cov(Y)). diff --git a/sourcecodes/bnt-master/KPMstats/README.txt b/sourcecodes/bnt-master/KPMstats/README.txt new file mode 100644 index 00000000..29e828d8 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/README.txt @@ -0,0 +1,4 @@ +KPMstats is a directory of miscellaneous statistics functions written by +Kevin Patrick Murphy and various other people (see individual file headers). + + diff --git a/sourcecodes/bnt-master/KPMstats/beta_sample.m b/sourcecodes/bnt-master/KPMstats/beta_sample.m new file mode 100644 index 00000000..82bacbe4 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/beta_sample.m @@ -0,0 +1,76 @@ +function r = betarnd(a,b,m,n); +%BETARND Random matrices from beta distribution. +% R = BETARND(A,B) returns a matrix of random numbers chosen +% from the beta distribution with parameters A and B. +% The size of R is the common size of A and B if both are matrices. +% If either parameter is a scalar, the size of R is the size of the other +% parameter. Alternatively, R = BETARND(A,B,M,N) returns an M by N matrix. + +% Reference: +% [1] L. Devroye, "Non-Uniform Random Variate Generation", +% Springer-Verlag, 1986 + +% Copyright (c) 1993-98 by The MathWorks, Inc. +% $Revision: 1.1.1.1 $ $Date: 2005/04/26 02:29:18 $ + +if nargin < 2, + error('Requires at least two input arguments'); +end + +if nargin == 2 + [errorcode rows columns] = rndcheck(2,2,a,b); +end + +if nargin == 3 + [errorcode rows columns] = rndcheck(3,2,a,b,m); +end + +if nargin == 4 + [errorcode rows columns] = rndcheck(4,2,a,b,m,n); +end + +if errorcode > 0 + error('Size information is inconsistent.'); +end + +r = zeros(rows,columns); + +% Use Theorem 4.1, case A (Devroye, page 430) to derive beta +% random numbers as a ratio of gamma random numbers. +if prod(size(a)) == 1 + a1 = a(ones(rows,1),ones(columns,1)); + g1 = gamrnd(a1,1); +else + g1 = gamrnd(a,1); +end +if prod(size(b)) == 1 + b1 = b(ones(rows,1),ones(columns,1)); + g2 = gamrnd(b1,1); +else + g2 = gamrnd(b,1); +end +r = g1 ./ (g1 + g2); + +% Return NaN if b is not positive. +if any(any(b <= 0)); + if prod(size(b) == 1) + tmp = NaN; + r = tmp(ones(rows,columns)); + else + k = find(b <= 0); + tmp = NaN; + r(k) = tmp(ones(size(k))); + end +end + +% Return NaN if a is not positive. +if any(any(a <= 0)); + if prod(size(a) == 1) + tmp = NaN; + r = tmp(ones(rows,columns)); + else + k = find(a <= 0); + tmp = NaN; + r(k) = tmp(ones(size(k))); + end +end diff --git a/sourcecodes/bnt-master/KPMstats/chisquared_histo.m b/sourcecodes/bnt-master/KPMstats/chisquared_histo.m new file mode 100644 index 00000000..1b649335 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/chisquared_histo.m @@ -0,0 +1,5 @@ +function s = chisquared_histo(h1, h2) +% Measure distance between 2 histograms (small numbers means more similar) +denom = h1 + h2; +denom = denom + (denom==0); +s = sum(((h1 - h2) .^ 2) ./ denom); diff --git a/sourcecodes/bnt-master/KPMstats/chisquared_prob.m b/sourcecodes/bnt-master/KPMstats/chisquared_prob.m new file mode 100644 index 00000000..4cf6fdc1 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/chisquared_prob.m @@ -0,0 +1,35 @@ +function P = chisquared_prob(X2,v) +%CHISQUARED_PROB computes the chi-squared probability function. +% P = CHISQUARED_PROB( X2, v ) returns P(X2|v), the probability +% of observing a chi-squared value <= X2 with v degrees of freedom. +% This is the probability that the sum of squares of v unit-variance +% normally-distributed random variables is <= X2. +% X2 and v may be matrices of the same size size, or either +% may be a scalar. +% +% e.g., CHISQUARED_PROB(5.99,2) returns 0.9500, verifying the +% 95% confidence bound for 2 degrees of freedom. This is also +% cross-checked in, e.g., Abramowitz & Stegun Table 26.8 +% +% See also CHISQUARED_TABLE +% +%Peter R. Shaw, WHOI + +% References: Press et al., Numerical Recipes, Cambridge, 1986; +% Abramowitz & Stegun, Handbook of Mathematical Functions, Dover, 1972. + +% Peter R. Shaw, Woods Hole Oceanographic Institution +% Woods Hole, MA 02543 +% (508) 457-2000 ext. 2473 pshaw@whoi.edu +% March, 1990; fixed Oct 1992 for version 4 + +% Computed using the Incomplete Gamma function, +% as given by Press et al. (Recipes) eq. (6.2.17) + +% Following nonsense is necessary from Matlab version 3 -> version 4 +versn_str=version; eval(['versn=' versn_str(1) ';']); +if versn<=3, %sigh + P = gamma(v/2, X2/2); +else + P = gammainc(X2/2, v/2); +end diff --git a/sourcecodes/bnt-master/KPMstats/chisquared_readme.txt b/sourcecodes/bnt-master/KPMstats/chisquared_readme.txt new file mode 100644 index 00000000..fc2ab64c --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/chisquared_readme.txt @@ -0,0 +1,36 @@ +From ftp://ftp.mathworks.com/pub/contrib/v4/stats/chisquared/ + +Chi-squared probability function. + +Abstract +m-files to compute the Chi-squared probability function, and +the percentage points of the probability function. + +P = CHISQUARED_PROB( X2, v ) returns P(X2|v), the probability +of observing a chi-squared value <= X2 with v degrees of freedom. +This is the probability that the sum of squares of v unit-variance +normally-distributed random variables is <= X2. + +Conversely: +X2 = CHISQUARED_TABLE( P, v ) returns the X2, the value of +chi-squared corresponding to v degrees of freedom and probability P. + +In reference textbooks, what is normally tabulated are the +percentage points of the chi-squared distribution; thus, one +would use CHISQUARED_TABLE rather than interpolate such a table. + +References: Press et al., Numerical Recipes, Cambridge, 1986; +Abramowitz & Stegun, Handbook of Mathematical Functions, +Dover, 1972; Table 26.8 + +Peter R. Shaw +Woods Hole Oceanographic Institution, Woods Hole, MA 02543 +(508) 457-2000 ext. 2473 +pshaw@whoi.edu + +--FILES GENERATED-- +README_chisq +chisquared_prob.m (function) -- chi-squared probability function +chisquared_table.m (function) -- "percentage points" (i.e., inverse) + of chi-squared probability function +chiaux.m (function) -- auxiliary fn. used by chisquared_table diff --git a/sourcecodes/bnt-master/KPMstats/chisquared_table.m b/sourcecodes/bnt-master/KPMstats/chisquared_table.m new file mode 100644 index 00000000..e9d191b3 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/chisquared_table.m @@ -0,0 +1,63 @@ +function X2 = chisquared_table(P,v) +%CHISQUARED_TABLE computes the "percentage points" of the +%chi-squared distribution, as in Abramowitz & Stegun Table 26.8 +% X2 = CHISQUARED_TABLE( P, v ) returns the value of chi-squared +% corresponding to v degrees of freedom and probability P. +% P is the probability that the sum of squares of v unit-variance +% normally-distributed random variables is <= X2. +% P and v may be matrices of the same size size, or either +% may be a scalar. +% +% e.g., to find the 95% confidence interval for 2 degrees +% of freedom, use CHISQUARED_TABLE( .95, 2 ), yielding 5.99, +% in agreement with Abramowitz & Stegun's Table 26.8 +% +% This result can be checked through the function +% CHISQUARED_PROB( 5.99, 2 ), yielding 0.9500 +% +% The familiar 1.96-sigma confidence bounds enclosing 95% of +% a 1-D gaussian is found through +% sqrt( CHISQUARED_TABLE( .95, 1 )), yielding 1.96 +% +% See also CHISQUARED_PROB +% +%Peter R. Shaw, WHOI +%Leslie Rosenfeld, MBARI + +% References: Press et al., Numerical Recipes, Cambridge, 1986; +% Abramowitz & Stegun, Handbook of Mathematical Functions, Dover, 1972. + +% Peter R. Shaw, Woods Hole Oceanographic Institution +% Woods Hole, MA 02543 pshaw@whoi.edu +% Leslie Rosenfeld, MBARI +% Last revision: Peter Shaw, Oct 1992: fsolve with version 4 + +% ** Calls function CHIAUX ** +% Computed using the Incomplete Gamma function, +% as given by Press et al. (Recipes) eq. (6.2.17) + +[mP,nP]=size(P); +[mv,nv]=size(v); +if mP~=mv | nP~=nv, + if mP==1 & nP==1, + P=P*ones(mv,nv); + elseif mv==1 & nv==1, + v=v*ones(mP,nP); + else + error('P and v must be the same size') + end +end +[m,n]=size(P); X2 = zeros(m,n); +for i=1:m, + for j=1:n, + if v(i,j)<=10, + x0=P(i,j)*v(i,j); + else + x0=v(i,j); + end +% Note: "old" and "new" calls to fsolve may or may not follow +% Matlab version 3.5 -> version 4 (so I'm keeping the old call around...) +% X2(i,j) = fsolve('chiaux',x0,zeros(16,1),[v(i,j),P(i,j)]); %(old call) + X2(i,j) = fsolve('chiaux',x0,zeros(16,1),[],[v(i,j),P(i,j)]); + end +end diff --git a/sourcecodes/bnt-master/KPMstats/clg_Mstep.m b/sourcecodes/bnt-master/KPMstats/clg_Mstep.m new file mode 100644 index 00000000..b7da903b --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/clg_Mstep.m @@ -0,0 +1,203 @@ +function [mu, Sigma, B] = clg_Mstep(w, Y, YY, YTY, X, XX, XY, varargin) +% MSTEP_CLG Compute ML/MAP estimates for a conditional linear Gaussian +% [mu, Sigma, B] = Mstep_clg(w, Y, YY, YTY, X, XX, XY, varargin) +% +% We fit P(Y|X,Q=i) = N(Y; B_i X + mu_i, Sigma_i) +% and w(i,t) = p(M(t)=i|y(t)) = posterior responsibility +% See www.ai.mit.edu/~murphyk/Papers/learncg.pdf. +% +% See process_options for how to specify the input arguments. +% +% INPUTS: +% w(i) = sum_t w(i,t) = responsibilities for each mixture component +% If there is only one mixture component (i.e., Q does not exist), +% then w(i) = N = nsamples, and +% all references to i can be replaced by 1. +% Y(:,i) = sum_t w(i,t) y(:,t) = weighted observations +% YY(:,:,i) = sum_t w(i,t) y(:,t) y(:,t)' = weighted outer product +% YTY(i) = sum_t w(i,t) y(:,t)' y(:,t) = weighted inner product +% You only need to pass in YTY if Sigma is to be estimated as spherical. +% +% In the regression context, we must also pass in the following +% X(:,i) = sum_t w(i,t) x(:,t) = weighted inputs +% XX(:,:,i) = sum_t w(i,t) x(:,t) x(:,t)' = weighted outer product +% XY(i) = sum_t w(i,t) x(:,t) y(:,t)' = weighted outer product +% +% Optional inputs (default values in []) +% +% 'cov_type' - 'full', 'diag' or 'spherical' ['full'] +% 'tied_cov' - 1 (Sigma) or 0 (Sigma_i) [0] +% 'clamped_cov' - pass in clamped value, or [] if unclamped [ [] ] +% 'clamped_mean' - pass in clamped value, or [] if unclamped [ [] ] +% 'clamped_weights' - pass in clamped value, or [] if unclamped [ [] ] +% 'cov_prior' - added to Sigma(:,:,i) to ensure psd [0.01*eye(d,d,Q)] +% +% If cov is tied, Sigma has size d*d. +% But diagonal and spherical covariances are represented in full size. + +[cov_type, tied_cov, ... + clamped_cov, clamped_mean, clamped_weights, cov_prior, ... + xs, ys, post] = ... + process_options(varargin, ... + 'cov_type', 'full', 'tied_cov', 0, 'clamped_cov', [], 'clamped_mean', [], ... + 'clamped_weights', [], 'cov_prior', [], ... + 'xs', [], 'ys', [], 'post', []); + +[Ysz Q] = size(Y); + +if isempty(X) % no regression + %B = []; + B2 = zeros(Ysz, 1, Q); + for i=1:Q + B(:,:,i) = B2(:,1:0,i); % make an empty array of size Ysz x 0 x Q + end + [mu, Sigma] = mixgauss_Mstep(w, Y, YY, YTY, varargin{:}); + return; +end + + +N = sum(w); +if isempty(cov_prior) + cov_prior = 0.01*repmat(eye(Ysz,Ysz), [1 1 Q]); +end +%YY = YY + cov_prior; % regularize the scatter matrix + +% Set any zero weights to one before dividing +% This is valid because w(i)=0 => Y(:,i)=0, etc +w = w + (w==0); + +Xsz = size(X,1); +% Append 1 to X to get Z +ZZ = zeros(Xsz+1, Xsz+1, Q); +ZY = zeros(Xsz+1, Ysz, Q); +for i=1:Q + ZZ(:,:,i) = [XX(:,:,i) X(:,i); + X(:,i)' w(i)]; + ZY(:,:,i) = [XY(:,:,i); + Y(:,i)']; +end + + +%%% Estimate mean and regression + +if ~isempty(clamped_weights) & ~isempty(clamped_mean) + B = clamped_weights; + mu = clamped_mean; +end +if ~isempty(clamped_weights) & isempty(clamped_mean) + B = clamped_weights; + % eqn 5 + mu = zeros(Ysz, Q); + for i=1:Q + mu(:,i) = (Y(:,i) - B(:,:,i)*X(:,i)) / w(i); + end +end +if isempty(clamped_weights) & ~isempty(clamped_mean) + mu = clamped_mean; + % eqn 3 + B = zeros(Ysz, Xsz, Q); + for i=1:Q + tmp = XY(:,:,i)' - mu(:,i)*X(:,i)'; + %B(:,:,i) = tmp * inv(XX(:,:,i)); + B(:,:,i) = (XX(:,:,i) \ tmp')'; + end +end +if isempty(clamped_weights) & isempty(clamped_mean) + mu = zeros(Ysz, Q); + B = zeros(Ysz, Xsz, Q); + % Nothing is clamped, so we must estimate B and mu jointly + for i=1:Q + % eqn 9 + if rcond(ZZ(:,:,i)) < 1e-10 + sprintf('clg_Mstep warning: ZZ(:,:,%d) is ill-conditioned', i); + % probably because there are too few cases for a high-dimensional input + ZZ(:,:,i) = ZZ(:,:,i) + 1e-5*eye(Xsz+1); + end + %A = ZY(:,:,i)' * inv(ZZ(:,:,i)); + A = (ZZ(:,:,i) \ ZY(:,:,i))'; + B(:,:,i) = A(:, 1:Xsz); + mu(:,i) = A(:, Xsz+1); + end +end + +if ~isempty(clamped_cov) + Sigma = clamped_cov; + return; +end + + +%%% Estimate covariance + +% Spherical +if cov_type(1)=='s' + if ~tied_cov + Sigma = zeros(Ysz, Ysz, Q); + for i=1:Q + % eqn 16 + A = [B(:,:,i) mu(:,i)]; + %s = trace(YTY(i) + A'*A*ZZ(:,:,i) - 2*A*ZY(:,:,i)) / (Ysz*w(i)); % wrong! + s = (YTY(i) + trace(A'*A*ZZ(:,:,i)) - trace(2*A*ZY(:,:,i))) / (Ysz*w(i)); + Sigma(:,:,i) = s*eye(Ysz,Ysz); + + %%%%%%%%%%%%%%%%%%% debug + if ~isempty(xs) + [nx T] = size(xs); + zs = [xs; ones(1,T)]; + yty = 0; + zAAz = 0; + yAz = 0; + for t=1:T + yty = yty + ys(:,t)'*ys(:,t) * post(i,t); + zAAz = zAAz + zs(:,t)'*A'*A*zs(:,t)*post(i,t); + yAz = yAz + ys(:,t)'*A*zs(:,t)*post(i,t); + end + assert(approxeq(yty, YTY(i))) + assert(approxeq(zAAz, trace(A'*A*ZZ(:,:,i)))) + assert(approxeq(yAz, trace(A*ZY(:,:,i)))) + s2 = (yty + zAAz - 2*yAz) / (Ysz*w(i)); + assert(approxeq(s,s2)) + end + %%%%%%%%%%%%%%% end debug + + end + else + S = 0; + for i=1:Q + % eqn 18 + A = [B(:,:,i) mu(:,i)]; + S = S + trace(YTY(i) + A'*A*ZZ(:,:,i) - 2*A*ZY(:,:,i)); + end + Sigma = repmat(S / (N*Ysz), [1 1 Q]); + end +else % Full/diagonal + if ~tied_cov + Sigma = zeros(Ysz, Ysz, Q); + for i=1:Q + A = [B(:,:,i) mu(:,i)]; + % eqn 10 + SS = (YY(:,:,i) - ZY(:,:,i)'*A' - A*ZY(:,:,i) + A*ZZ(:,:,i)*A') / w(i); + if cov_type(1)=='d' + Sigma(:,:,i) = diag(diag(SS)); + else + Sigma(:,:,i) = SS; + end + end + else % tied + SS = zeros(Ysz, Ysz); + for i=1:Q + A = [B(:,:,i) mu(:,i)]; + % eqn 13 + SS = SS + (YY(:,:,i) - ZY(:,:,i)'*A' - A*ZY(:,:,i) + A*ZZ(:,:,i)*A'); + end + SS = SS / N; + if cov_type(1)=='d' + Sigma = diag(diag(SS)); + else + Sigma = SS; + end + Sigma = repmat(Sigma, [1 1 Q]); + end +end + +Sigma = Sigma + cov_prior; + diff --git a/sourcecodes/bnt-master/KPMstats/clg_Mstep_simple.m b/sourcecodes/bnt-master/KPMstats/clg_Mstep_simple.m new file mode 100644 index 00000000..5b66bb0e --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/clg_Mstep_simple.m @@ -0,0 +1,52 @@ +function [mu, B] = clg_Mstep_simple(w, Y, YY, YTY, X, XX, XY) +% CLG_MSTEP_SIMPLE Same as CLG_MSTEP, but doesn;t estimate Sigma, so is slightly faster +% function [mu, B] = clg_Mstep_simple(w, Y, YY, YTY, X, XX, XY) +% +% See clg_Mstep for details. +% Unlike clg_Mstep, there are no optional arguments, which are slow to process +% if this function is inside a tight loop. + +[Ysz Q] = size(Y); + +if isempty(X) % no regression + %B = []; + B2 = zeros(Ysz, 1, Q); + for i=1:Q + B(:,:,i) = B2(:,1:0,i); % make an empty array of size Ysz x 0 x Q + end + [mu, Sigma] = mixgauss_Mstep(w, Y, YY, YTY); + return; +end + +N = sum(w); +%YY = YY + cov_prior; % regularize the scatter matrix + +% Set any zero weights to one before dividing +% This is valid because w(i)=0 => Y(:,i)=0, etc +w = w + (w==0); + +Xsz = size(X,1); +% Append 1 to X to get Z +ZZ = zeros(Xsz+1, Xsz+1, Q); +ZY = zeros(Xsz+1, Ysz, Q); +for i=1:Q + ZZ(:,:,i) = [XX(:,:,i) X(:,i); + X(:,i)' w(i)]; + ZY(:,:,i) = [XY(:,:,i); + Y(:,i)']; +end + +mu = zeros(Ysz, Q); +B = zeros(Ysz, Xsz, Q); +for i=1:Q + % eqn 9 + if rcond(ZZ(:,:,i)) < 1e-10 + sprintf('clg_Mstep warning: ZZ(:,:,%d) is ill-conditioned', i); + %probably because there are too few cases for a high-dimensional input + ZZ(:,:,i) = ZZ(:,:,i) + 1e-5*eye(Xsz+1); + end + %A = ZY(:,:,i)' * inv(ZZ(:,:,i)); + A = (ZZ(:,:,i) \ ZY(:,:,i))'; + B(:,:,i) = A(:, 1:Xsz); + mu(:,i) = A(:, Xsz+1); +end diff --git a/sourcecodes/bnt-master/KPMstats/clg_prob.m b/sourcecodes/bnt-master/KPMstats/clg_prob.m new file mode 100644 index 00000000..3ad7b761 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/clg_prob.m @@ -0,0 +1,14 @@ +function p = eval_pdf_clg(X,Y,mu,Sigma,W) +% function p = eval_pdf_clg(X,Y,mu,Sigma,W) +% +% p(c,t) = N(Y(:,t); mu(:,c) + W(:,:,c)*X(:,t), Sigma(:,:,c)) + +[d T] = size(Y); +[d nc] = size(mu); +p = zeros(nc,T); +for c=1:nc + denom = (2*pi)^(d/2)*sqrt(abs(det(Sigma(:,:,c)))); + M = repmat(mu(:,c), 1, T) + W(:,:,c)*X; + mahal = sum(((Y-M)'*inv(Sigma(:,:,c))).*(Y-M)',2); + p(c,:) = (exp(-0.5*mahal) / denom)'; +end diff --git a/sourcecodes/bnt-master/KPMstats/condGaussToJoint.m b/sourcecodes/bnt-master/KPMstats/condGaussToJoint.m new file mode 100644 index 00000000..59075675 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/condGaussToJoint.m @@ -0,0 +1,22 @@ +function [muXY, SigmaXY] = condGaussToJoint(muX, SigmaX, muY, SigmaY, WYgivenX) + +% Compute P(X,Y) from P(X) * P(Y|X) where P(X)=N(X;muX,SigmaX) +% and P(Y|X) = N(Y; WX + muY, SigmaY) + +% For details on how to compute a Gaussian from a Bayes net +% - "Gaussian Influence Diagrams", R. Shachter and C. R. Kenley, Management Science, 35(5):527--550, 1989. + +% size(W) = dy x dx +dx = length(muX); +dy = length(muY); +muXY = [muX(:); WYgivenX*muX(:) + muY]; + +W = [zeros(dx,dx) WYgivenX'; + zeros(dy,dx) zeros(dy,dy)]; +D = [SigmaX zeros(dx,dy); + zeros(dy,dx) SigmaY]; + +U = inv(eye(size(W)) - W')'; +SigmaXY = U' * D * U; + + diff --git a/sourcecodes/bnt-master/KPMstats/cond_indep_fisher_z.m b/sourcecodes/bnt-master/KPMstats/cond_indep_fisher_z.m new file mode 100644 index 00000000..e116ba59 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/cond_indep_fisher_z.m @@ -0,0 +1,142 @@ +function [CI, r, p] = cond_indep_fisher_z(X, Y, S, C, N, alpha) +% COND_INDEP_FISHER_Z Test if X indep Y given Z using Fisher's Z test +% CI = cond_indep_fisher_z(X, Y, S, C, N, alpha) +% +% C is the covariance (or correlation) matrix +% N is the sample size +% alpha is the significance level (default: 0.05) +% +% See p133 of T. Anderson, "An Intro. to Multivariate Statistical Analysis", 1984 + +if nargin < 6, alpha = 0.05; end + +r = partial_corr_coef(C, X, Y, S); +z = 0.5*log( (1+r)/(1-r) ); +z0 = 0; +W = sqrt(N - length(S) - 3)*(z-z0); % W ~ N(0,1) +cutoff = norminv(1 - 0.5*alpha); % P(|W| <= cutoff) = 0.95 +%cutoff = mynorminv(1 - 0.5*alpha); % P(|W| <= cutoff) = 0.95 +if abs(W) < cutoff + CI = 1; +else % reject the null hypothesis that rho = 0 + CI = 0; +end +p = normcdf(W); +%p = mynormcdf(W); + +%%%%%%%%% + +function p = normcdf(x,mu,sigma) +%NORMCDF Normal cumulative distribution function (cdf). +% P = NORMCDF(X,MU,SIGMA) computes the normal cdf with mean MU and +% standard deviation SIGMA at the values in X. +% +% The size of P is the common size of X, MU and SIGMA. A scalar input +% functions as a constant matrix of the same size as the other inputs. +% +% Default values for MU and SIGMA are 0 and 1 respectively. + +% References: +% [1] M. Abramowitz and I. A. Stegun, "Handbook of Mathematical +% Functions", Government Printing Office, 1964, 26.2. + +% Copyright (c) 1993-98 by The MathWorks, Inc. +% $Revision: 1.1.1.1 $ $Date: 2005/04/26 02:29:18 $ + +if nargin < 3, + sigma = 1; +end + +if nargin < 2; + mu = 0; +end + +[errorcode x mu sigma] = distchck(3,x,mu,sigma); + +if errorcode > 0 + error('Requires non-scalar arguments to match in size.'); +end + +% Initialize P to zero. +p = zeros(size(x)); + +% Return NaN if SIGMA is not positive. +k1 = find(sigma <= 0); +if any(k1) + tmp = NaN; + p(k1) = tmp(ones(size(k1))); +end + +% Express normal CDF in terms of the error function. +k = find(sigma > 0); +if any(k) + p(k) = 0.5 * erfc( - (x(k) - mu(k)) ./ (sigma(k) * sqrt(2))); +end + +% Make sure that round-off errors never make P greater than 1. +k2 = find(p > 1); +if any(k2) + p(k2) = ones(size(k2)); +end + +%%%%%%%% + +function x = norminv(p,mu,sigma); +%NORMINV Inverse of the normal cumulative distribution function (cdf). +% X = NORMINV(P,MU,SIGMA) finds the inverse of the normal cdf with +% mean, MU, and standard deviation, SIGMA. +% +% The size of X is the common size of the input arguments. A scalar input +% functions as a constant matrix of the same size as the other inputs. +% +% Default values for MU and SIGMA are 0 and 1 respectively. + +% References: +% [1] M. Abramowitz and I. A. Stegun, "Handbook of Mathematical +% Functions", Government Printing Office, 1964, 7.1.1 and 26.2.2 + +% Copyright (c) 1993-98 by The MathWorks, Inc. +% $Revision: 1.1.1.1 $ $Date: 2005/04/26 02:29:18 $ + +if nargin < 3, + sigma = 1; +end + +if nargin < 2; + mu = 0; +end + +[errorcode p mu sigma] = distchck(3,p,mu,sigma); + +if errorcode > 0 + error('Requires non-scalar arguments to match in size.'); +end + +% Allocate space for x. +x = zeros(size(p)); + +% Return NaN if the arguments are outside their respective limits. +k = find(sigma <= 0 | p < 0 | p > 1); +if any(k) + tmp = NaN; + x(k) = tmp(ones(size(k))); +end + +% Put in the correct values when P is either 0 or 1. +k = find(p == 0); +if any(k) + tmp = Inf; + x(k) = -tmp(ones(size(k))); +end + +k = find(p == 1); +if any(k) + tmp = Inf; + x(k) = tmp(ones(size(k))); +end + +% Compute the inverse function for the intermediate values. +k = find(p > 0 & p < 1 & sigma > 0); +if any(k), + x(k) = sqrt(2) * sigma(k) .* erfinv(2 * p(k) - 1) + mu(k); +end diff --git a/sourcecodes/bnt-master/KPMstats/condgaussTrainObserved.m b/sourcecodes/bnt-master/KPMstats/condgaussTrainObserved.m new file mode 100644 index 00000000..2dd4cafa --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/condgaussTrainObserved.m @@ -0,0 +1,27 @@ +function [mu, Sigma] = mixgaussTrainObserved(obsData, hiddenData, nstates, varargin); +% mixgaussTrainObserved Max likelihood estimates of conditional Gaussian from raw data +% function [mu, Sigma] = mixgaussTrainObserved(obsData, hiddenData, nstates, ...); +% +% Input: +% obsData(:,i) +% hiddenData(i) - this is the mixture component label for example i +% Optional arguments - same as mixgauss_Mstep +% +% Output: +% mu(:,q) +% Sigma(:,:,q) - same as mixgauss_Mstep + +[D numex] = size(obsData); +Y = zeros(D, nstates); +YY = zeros(D,D,nstates); +YTY = zeros(nstates,1); +w = zeros(nstates, 1); +for q=1:nstates + ndx = find(hiddenData==q); + w(q) = length(ndx); % each data point has probability 1 of being in this cluster + data = obsData(:,ndx); + Y(:,q) = sum(data,2); + YY(:,:,q) = data*data'; + YTY(q) = sum(diag(data'*data)); +end +[mu, Sigma] = mixgauss_Mstep(w, Y, YY, YTY, varargin{:}); diff --git a/sourcecodes/bnt-master/KPMstats/condgauss_sample.m b/sourcecodes/bnt-master/KPMstats/condgauss_sample.m new file mode 100644 index 00000000..17ca0eff --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/condgauss_sample.m @@ -0,0 +1,11 @@ +function x = mixgauss_sample(mu, Sigma, labels) +% MIXGAUSS_SAMPLE Sample from a mixture of Gaussians given known mixture labels +% function x = mixgauss_sample(mu, Sigma, labels) + +T = length(labels); +[D Q] = size(mu); +x = zeros(D,T); +for q=1:Q + ndx = find(labels==q); + x(:,ndx) = gaussian_sample(mu(:,q)', Sigma(:,:,q), length(ndx))'; +end diff --git a/sourcecodes/bnt-master/KPMstats/convertBinaryLabels.m b/sourcecodes/bnt-master/KPMstats/convertBinaryLabels.m new file mode 100644 index 00000000..aeb1e1ee --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/convertBinaryLabels.m @@ -0,0 +1,3 @@ +% labels01 = (labelsPM+1)/2; % maps -1->0, +1->1 +% labelsPM = (2*labels01)-1; % maps 0,1 -> -1,1 + diff --git a/sourcecodes/bnt-master/KPMstats/cwr_demo.m b/sourcecodes/bnt-master/KPMstats/cwr_demo.m new file mode 100644 index 00000000..75090535 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/cwr_demo.m @@ -0,0 +1,124 @@ +% Compare my code with +% http://www.media.mit.edu/physics/publications/books/nmm/files/index.html +% +% cwm.m +% (c) Neil Gershenfeld 9/1/97 +% 1D Cluster-Weighted Modeling example +% +clear all +figure; +seed = 0; +rand('state', seed); +randn('state', seed); +x = (-10:10)'; +y = (x > 0); +npts = length(x); +plot(x,y,'+') +xlabel('x') +ylabel('y') +nclusters = 4; +nplot = 100; +xplot = 24*(1:nplot)'/nplot - 12; + +mux = 20*rand(1,nclusters) - 10; +muy = zeros(1,nclusters); +varx = ones(1,nclusters); +vary = ones(1,nclusters); +pc = 1/nclusters * ones(1,nclusters); +niterations = 5; +eps = 0.01; + + +I = repmat(eye(1,1), [1 1 nclusters]); +O = repmat(zeros(1,1), [1 1 nclusters]); +X = x(:)'; +Y = y(:)'; + +cwr = cwr_em(X, Y, nclusters, 'muX', mux, 'muY', muy, 'SigmaX', I, ... + 'cov_typeX', 'spherical', 'SigmaY', I, 'cov_typeY', 'spherical', ... + 'priorC', pc, 'weightsY', O, 'create_init_params', 0, ... + 'clamp_weights', 1, 'max_iter', niterations, ... + 'cov_priorX', eps*ones(1,1,nclusters), ... + 'cov_priorY', eps*ones(1,1,nclusters)); + + +% Gershenfeld's EM code +for step = 1:niterations + pplot = exp(-(kron(xplot,ones(1,nclusters)) ... + - kron(ones(nplot,1),mux)).^2 ... + ./ (2*kron(ones(nplot,1),varx))) ... + ./ sqrt(2*pi*kron(ones(nplot,1),varx)) ... + .* kron(ones(nplot,1),pc); + plot(xplot,pplot,'k'); + pause(0); + px = exp(-(kron(x,ones(1,nclusters)) ... + - kron(ones(npts,1),mux)).^2 ... + ./ (2*kron(ones(npts,1),varx))) ... + ./ sqrt(2*pi*kron(ones(npts,1),varx)); + py = exp(-(kron(y,ones(1,nclusters)) ... + - kron(ones(npts,1),muy)).^2 ... + ./ (2*kron(ones(npts,1),vary))) ... + ./ sqrt(2*pi*kron(ones(npts,1),vary)); + p = px .* py .* kron(ones(npts,1),pc); + pp = p ./ kron(sum(p,2),ones(1,nclusters)); + pc = sum(pp)/npts; + yfit = sum(kron(ones(npts,1),muy) .* p,2) ... + ./ sum(p,2); + mux = sum(kron(x,ones(1,nclusters)) .* pp) ... + ./ (npts*pc); + varx = eps + sum((kron(x,ones(1,nclusters)) ... + - kron(ones(npts,1),mux)).^2 .* pp) ... + ./ (npts*pc); + muy = sum(kron(y,ones(1,nclusters)) .* pp) ... + ./ (npts*pc); + vary = eps + sum((kron(y,ones(1,nclusters)) ... + - kron(ones(npts,1),muy)).^2 .* pp) ... + ./ (npts*pc); +end + + +% Check equal +cwr_pc = cwr.priorC'; +assert(approxeq(cwr_pc, pc)) +cwr_mux = cwr.muX; +assert(approxeq(mux, cwr_mux)) +cwr_SigmaX = squeeze(cwr.SigmaX)'; +assert(approxeq(varx, cwr_SigmaX)) +cwr_muy = cwr.muY; +assert(approxeq(muy, cwr_muy)) +cwr_SigmaY = squeeze(cwr.SigmaY)'; +assert(approxeq(vary, cwr_SigmaY)) + + +% Prediction + +X = xplot(:)'; +[cwr_mu, Sigma, post] = cwr_predict(cwr, X); +cwr_ystd = squeeze(Sigma)'; + +% pplot(t,c) +pplot = exp(-(kron(xplot,ones(1,nclusters)) ... + - kron(ones(nplot,1),mux)).^2 ... + ./ (2*kron(ones(nplot,1),varx))) ... + ./ sqrt(2*pi*kron(ones(nplot,1),varx)) ... + .* kron(ones(nplot,1),pc); +yplot = sum(kron(ones(nplot,1),muy) .* pplot,2) ... + ./ sum(pplot,2); +ystdplot = sum(kron(ones(nplot,1),(muy.^2+vary)) .* pplot,2) ... + ./ sum(pplot,2) - yplot.^2; + + +% Check equal +assert(approxeq(yplot(:)', cwr_mu(:)')) +assert(approxeq(ystdplot, cwr_ystd)) +assert(approxeq(pplot ./ repmat(sum(pplot,2), 1, nclusters),post') ) + +plot(xplot,yplot,'k'); +hold on +plot(xplot,yplot+ystdplot,'k--'); +plot(xplot,yplot-ystdplot,'k--'); +plot(x,y,'k+'); +axis([-12 12 -1 1.1]); +plot(xplot,.8*pplot/max(max(pplot))-1,'k') +hold off + diff --git a/sourcecodes/bnt-master/KPMstats/cwr_em.m b/sourcecodes/bnt-master/KPMstats/cwr_em.m new file mode 100644 index 00000000..96bf491e --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/cwr_em.m @@ -0,0 +1,161 @@ +function cwr = cwr_em(X, Y, nc, varargin) +% CWR_LEARN Fit the parameters of a cluster weighted regression model using EM +% function cwr = cwr_learn(X, Y, ...) +% +% X(:, t) is the t'th input example +% Y(:, t) is the t'th output example +% nc is the number of clusters +% +% Kevin Murphy, May 2003 + +[max_iter, thresh, cov_typeX, cov_typeY, clamp_weights, ... + muX, muY, SigmaX, SigmaY, weightsY, priorC, create_init_params, ... +cov_priorX, cov_priorY, verbose, regress, clamp_covX, clamp_covY] = process_options(... + varargin, 'max_iter', 10, 'thresh', 1e-2, 'cov_typeX', 'full', ... + 'cov_typeY', 'full', 'clamp_weights', 0, ... + 'muX', [], 'muY', [], 'SigmaX', [], 'SigmaY', [], 'weightsY', [], 'priorC', [], ... + 'create_init_params', 1, 'cov_priorX', [], 'cov_priorY', [], 'verbose', 0, ... + 'regress', 1, 'clamp_covX', 0, 'clamp_covY', 0); + +[nx N] = size(X); +[ny N2] = size(Y); +if N ~= N2 + error(sprintf('nsamples X (%d) ~= nsamples Y (%d)', N, N2)); +end +%if N < nx +% fprintf('cwr_em warning: dim X (%d) > nsamples X (%d)\n', nx, N); +%end +if (N < nx) & regress + fprintf('cwr_em warning: dim X = %d, nsamples X = %d\n', nx, N); +end +if (N < ny) + fprintf('cwr_em warning: dim Y = %d, nsamples Y = %d\n', ny, N); +end +if (nc > N) + error(sprintf('cwr_em: more centers (%d) than data', nc)) +end + +if nc==1 + % No latent variable, so there is a closed-form solution + w = 1/N; + WYbig = Y*w; + WYY = WYbig * Y'; + WY = sum(WYbig, 2); + WYTY = sum(diag(WYbig' * Y)); + cwr.priorC = 1; + cwr.SigmaX = []; + if ~regress + % This is just fitting an unconditional Gaussian + cwr.weightsY = []; + [cwr.muY, cwr.SigmaY] = ... + mixgauss_Mstep(1, WY, WYY, WYTY, ... + 'cov_type', cov_typeY, 'cov_prior', cov_priorY); + % There is a much easier way... + assert(approxeq(cwr.muY, mean(Y'))) + assert(approxeq(cwr.SigmaY, cov(Y') + 0.01*eye(ny))) + else + % This is just linear regression + WXbig = X*w; + WXX = WXbig * X'; + WX = sum(WXbig, 2); + WXTX = sum(diag(WXbig' * X)); + WXY = WXbig * Y'; + [cwr.muY, cwr.SigmaY, cwr.weightsY] = ... + clg_Mstep(1, WY, WYY, WYTY, WX, WXX, WXY, ... + 'cov_type', cov_typeY, 'cov_prior', cov_priorY); + end + if clamp_covY, cwr.SigmaY = SigmaY; end + if clamp_weights, cwr.weightsY = weightsY; end + return; +end + + +if create_init_params + [cwr.muX, cwr.SigmaX] = mixgauss_init(nc, X, cov_typeX); + [cwr.muY, cwr.SigmaY] = mixgauss_init(nc, Y, cov_typeY); + cwr.weightsY = zeros(ny, nx, nc); + cwr.priorC = normalize(ones(nc,1)); +else + cwr.muX = muX; cwr.muY = muY; cwr.SigmaX = SigmaX; cwr.SigmaY = SigmaY; + cwr.weightsY = weightsY; cwr.priorC = priorC; +end + + +if clamp_covY, cwr.SigmaY = SigmaY; end +if clamp_covX, cwr.SigmaX = SigmaX; end +if clamp_weights, cwr.weightsY = weightsY; end + +previous_loglik = -inf; +num_iter = 1; +converged = 0; + +while (num_iter <= max_iter) & ~converged + + % E step + + [likXandY, likYgivenX, post] = cwr_prob(cwr, X, Y); + loglik = sum(log(likXandY)); + % extract expected sufficient statistics + w = sum(post,2); % post(c,t) + WYY = zeros(ny, ny, nc); + WY = zeros(ny, nc); + WYTY = zeros(nc,1); + + WXX = zeros(nx, nx, nc); + WX = zeros(nx, nc); + WXTX = zeros(nc, 1); + WXY = zeros(nx,ny,nc); + %WYY = repmat(reshape(w, [1 1 nc]), [ny ny 1]) .* repmat(Y*Y', [1 1 nc]); + for c=1:nc + weights = repmat(post(c,:), ny, 1); + WYbig = Y .* weights; + WYY(:,:,c) = WYbig * Y'; + WY(:,c) = sum(WYbig, 2); + WYTY(c) = sum(diag(WYbig' * Y)); + + weights = repmat(post(c,:), nx, 1); % weights(nx, nsamples) + WXbig = X .* weights; + WXX(:,:,c) = WXbig * X'; + WX(:,c) = sum(WXbig, 2); + WXTX(c) = sum(diag(WXbig' * X)); + WXY(:,:,c) = WXbig * Y'; + end + + % M step + % Q -> X is called Q->Y in Mstep_clg + [cwr.muX, cwr.SigmaX] = mixgauss_Mstep(w, WX, WXX, WXTX, ... + 'cov_type', cov_typeX, 'cov_prior', cov_priorX); + for c=1:nc + assert(is_psd(cwr.SigmaX(:,:,c))) + end + + if clamp_weights % affects estimate of mu and Sigma + W = cwr.weightsY; + else + W = []; + end + [cwr.muY, cwr.SigmaY, cwr.weightsY] = ... + clg_Mstep(w, WY, WYY, WYTY, WX, WXX, WXY, ... + 'cov_type', cov_typeY, 'clamped_weights', W, ... + 'cov_prior', cov_priorY); + %'xs', X, 'ys', Y, 'post', post); % debug + %a = linspace(min(Y(2,:)), max(Y(2,:)), nc+2); + %cwr.muY(2,:) = a(2:end-1); + + cwr.priorC = normalize(w); + + for c=1:nc + assert(is_psd(cwr.SigmaY(:,:,c))) + end + + if clamp_covY, cwr.SigmaY = SigmaY; end + if clamp_covX, cwr.SigmaX = SigmaX; end + if clamp_weights, cwr.weightsY = weightsY; end + + if verbose, fprintf(1, 'iteration %d, loglik = %f\n', num_iter, loglik); end + num_iter = num_iter + 1; + converged = em_converged(loglik, previous_loglik, thresh); + previous_loglik = loglik; + +end + diff --git a/sourcecodes/bnt-master/KPMstats/cwr_predict.m b/sourcecodes/bnt-master/KPMstats/cwr_predict.m new file mode 100644 index 00000000..9fe54cab --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/cwr_predict.m @@ -0,0 +1,57 @@ +function [mu, Sigma, weights, mask] = cwr_predict(cwr, X, mask_data) +% CWR_PREDICT cluster weighted regression: predict Y given X +% function [mu, Sigma] = cwr_predict(cwr, X) +% +% mu(:,t) = E[Y|X(:,t)] = sum_c P(c | X(:,t)) E[Y|c, X(:,t)] +% Sigma(:,:,t) = Cov[Y|X(:,t)] +% +% [mu, Sigma, weights, mask] = cwr_predict(cwr, X, mask_data) +% mask(i) = sum_t sum_c p(mask_data(:,i) | X(:,t), c) P(c|X(:,t)) +% This evaluates the predictive density on a set of points +% (This is only sensible if T=1, ie. X is a single vector) + +[nx T] = size(X); +[ny nx nc] = size(cwr.weightsY); +mu = zeros(ny, T); +Sigma = zeros(ny, ny, T); + +if nargout == 4 + comp_mask = 1; + N = size(mask_data,2); + mask = zeros(N,1); +else + comp_mask = 0; +end + +if nc==1 + if isempty(cwr.weightsY) + mu = repmat(cwr.muY, 1, T); + Sigma = repmat(cwr.SigmaY, [1 1 T]); + else + mu = repmat(cwr.muY, 1, T) + cwr.weightsY * X; + Sigma = repmat(cwr.SigmaY, [1 1 T]); + %for t=1:T + % mu(:,t) = cwr.muY + cwr.weightsY*X(:,t); + % Sigma(:,:,t) = cwr.SigmaY; + %end + end + if comp_mask, mask = gaussian_prob(mask_data, mu, Sigma); end + weights = []; + return; +end + + +% likX(c,t) = p(x(:,t) | c) +likX = mixgauss_prob(X, cwr.muX, cwr.SigmaX); +weights = normalize(repmat(cwr.priorC, 1, T) .* likX, 1); +for t=1:T + mut = zeros(ny, nc); + for c=1:nc + mut(:,c) = cwr.muY(:,c) + cwr.weightsY(:,:,c)*X(:,t); + if comp_mask + mask = mask + gaussian_prob(mask_data, mut(:,c), cwr.SigmaY(:,:,c)) * weights(c); + end + end + %w = normalise(cwr.priorC(:) .* likX(:,t)); + [mu(:,t), Sigma(:,:,t)] = collapse_mog(mut, cwr.SigmaY, weights(:,t)); +end diff --git a/sourcecodes/bnt-master/KPMstats/cwr_prob.m b/sourcecodes/bnt-master/KPMstats/cwr_prob.m new file mode 100644 index 00000000..6971e0de --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/cwr_prob.m @@ -0,0 +1,39 @@ +function [likXandY, likYgivenX, post] = cwr_prob(cwr, X, Y); +% CWR_EVAL_PDF cluster weighted regression: evaluate likelihood of Y given X +% function [likXandY, likYgivenX, post] = cwr_prob(cwr, X, Y); +% +% likXandY(t) = p(x(:,t), y(:,t)) +% likXgivenY(t) = p(x(:,t)| y(:,t)) +% post(c,t) = p(c | x(:,t), y(:,t)) + +[nx N] = size(X); +nc = length(cwr.priorC); + +if nc == 1 + [mu, Sigma] = cwr_predict(cwr, X); + likY = gaussian_prob(Y, mu, Sigma); + likXandY = likY; + likYgivenX = likY; + post = ones(1,N); + return; +end + + +% likY(c,t) = p(y(:,t) | c) +likY = clg_prob(X, Y, cwr.muY, cwr.SigmaY, cwr.weightsY); + +% likX(c,t) = p(x(:,t) | c) +[junk, likX] = mixgauss_prob(X, cwr.muX, cwr.SigmaX); +likX = squeeze(likX); + +% prior(c,t) = p(c) +prior = repmat(cwr.priorC(:), 1, N); + +post = likX .* likY .* prior; +likXandY = sum(post, 1); +post = post ./ repmat(likXandY, nc, 1); +%loglik = sum(log(lik)); +%loglik = log(lik); + +likX = sum(likX .* prior, 1); +likYgivenX = likXandY ./ likX; diff --git a/sourcecodes/bnt-master/KPMstats/cwr_readme.txt b/sourcecodes/bnt-master/KPMstats/cwr_readme.txt new file mode 100644 index 00000000..28f6ce72 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/cwr_readme.txt @@ -0,0 +1,20 @@ +This directory implements Cluster Weighted Regression, as described in +Neil Gershenfeld, "The nature of mathematical modelling", p182. +(See also http://www.media.mit.edu/physics/publications/books/nmm/files/index.html) + +Written by K. Murphy, 2 May 2003 + +The model is as follows: + +X<--| +| Q +v | +Y<-- + +where Q is a discrete latent mixture variable. + +A mixture of experts has an X->Q arc instead of a Q->X arc; +the X->Q arc is modelled by a softmax, which is slightly harder to fit than a +mixture of Gaussians. + + diff --git a/sourcecodes/bnt-master/KPMstats/cwr_test.m b/sourcecodes/bnt-master/KPMstats/cwr_test.m new file mode 100644 index 00000000..dce1ed65 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/cwr_test.m @@ -0,0 +1,80 @@ +% Verify that my code gives the same results as the 1D example at +% http://www.media.mit.edu/physics/publications/books/nmm/files/cwm.m + +seed = 0; +rand('state', seed); +randn('state', seed); +x = (-10:10)'; +y = double(x > 0); +npts = length(x); +plot(x,y,'+') + +nclusters = 4; +nplot = 100; +xplot = 24*(1:nplot)'/nplot - 12; + +mux = 20*rand(1,nclusters) - 10; +muy = zeros(1,nclusters); +varx = ones(1,nclusters); +vary = ones(1,nclusters); +pc = 1/nclusters * ones(1,nclusters); + + +I = repmat(eye(1,1), [1 1 nclusters]); +O = repmat(zeros(1,1), [1 1 nclusters]); +X = x(:)'; +Y = y(:)'; + +% Do 1 iteration of EM + +%cwr = cwr_em(X, Y, nclusters, 'muX', mux, 'muY', muy, 'SigmaX', I, 'cov_typeX', 'spherical', 'SigmaY', I, 'cov_typeY', 'spherical', 'priorC', pc, 'weightsY', O, 'init_params', 0, 'clamp_weights', 1, 'max_iter', 1, 'cov_priorX', zeros(1,1,nclusters), 'cov_priorY', zeros(1,1,nclusters)); + +cwr = cwr_em(X, Y, nclusters, 'muX', mux, 'muY', muy, 'SigmaX', I, 'cov_typeX', 'spherical', 'SigmaY', I, 'cov_typeY', 'spherical', 'priorC', pc, 'weightsY', O, 'create_init_params', 0, 'clamp_weights', 1, 'max_iter', 1); + + +% Check this matches Gershenfeld's code + +% E step +% px(t,c) = prob(x(t) | c) +px = exp(-(kron(x,ones(1,nclusters)) ... + - kron(ones(npts,1),mux)).^2 ... + ./ (2*kron(ones(npts,1),varx))) ... + ./ sqrt(2*pi*kron(ones(npts,1),varx)); +py = exp(-(kron(y,ones(1,nclusters)) ... + - kron(ones(npts,1),muy)).^2 ... + ./ (2*kron(ones(npts,1),vary))) ... + ./ sqrt(2*pi*kron(ones(npts,1),vary)); +p = px .* py .* kron(ones(npts,1),pc); +pp = p ./ kron(sum(p,2),ones(1,nclusters)); + +% M step +eps = 0.01; +pc2 = sum(pp)/npts; + +mux2 = sum(kron(x,ones(1,nclusters)) .* pp) ... + ./ (npts*pc2); +varx2 = eps + sum((kron(x,ones(1,nclusters)) ... + - kron(ones(npts,1),mux2)).^2 .* pp) ... + ./ (npts*pc2); +muy2 = sum(kron(y,ones(1,nclusters)) .* pp) ... + ./ (npts*pc2); +vary2 = eps + sum((kron(y,ones(1,nclusters)) ... + - kron(ones(npts,1),muy2)).^2 .* pp) ... + ./ (npts*pc2); + + +denom = (npts*pc2); +% denom(c) = N*pc(c) = w(c) = sum_t pp(c,t) +% since pc(c) = sum_t pp(c,t) / N + +cwr_mux = cwr.muX; +assert(approxeq(mux2, cwr_mux)) +cwr_SigmaX = squeeze(cwr.SigmaX)'; +assert(approxeq(varx2, cwr_SigmaX)) + +cwr_muy = cwr.muY; +assert(approxeq(muy2, cwr_muy)) +cwr_SigmaY = squeeze(cwr.SigmaY)'; +assert(approxeq(vary2, cwr_SigmaY)) + + diff --git a/sourcecodes/bnt-master/KPMstats/dirichlet_sample.m b/sourcecodes/bnt-master/KPMstats/dirichlet_sample.m new file mode 100644 index 00000000..2ab09afb --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/dirichlet_sample.m @@ -0,0 +1,18 @@ +function theta = dirichlet_sample(alpha, N) +% SAMPLE_DIRICHLET Sample N vectors from Dir(alpha(1), ..., alpha(k)) +% theta = sample_dirichlet(alpha, N) +% theta(i,j) = i'th sample of theta_j, where theta ~ Dir + +% We use the method from p. 482 of "Bayesian Data Analysis", Gelman et al. + +assert(alpha > 0); +k = length(alpha); +theta = zeros(N, k); +scale = 1; % arbitrary +for i=1:k + %theta(:,i) = gamrnd(alpha(i), scale, N, 1); + theta(:,i) = gamma_sample(alpha(i), scale, N, 1); +end +%theta = mk_stochastic(theta); +S = sum(theta,2); +theta = theta ./ repmat(S, 1, k); diff --git a/sourcecodes/bnt-master/KPMstats/dirichletpdf.m b/sourcecodes/bnt-master/KPMstats/dirichletpdf.m new file mode 100644 index 00000000..8b789844 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/dirichletpdf.m @@ -0,0 +1,39 @@ +function p = dirichletpdf(x, alpha) +%DIRICHLETPDF Dirichlet probability density function. +% p = dirichletpdf(x, alpha) returns the probability of vector +% x under the Dirichlet distribution with parameter vector +% alpha. +% +% Author: David Ross + +%------------------------------------------------- +% Check the input +%------------------------------------------------- +error(nargchk(2,2,nargin)); + +% enusre alpha is a vector +if min(size(alpha)) ~= 1 | ndims(alpha) > 2 | length(alpha) == 1 + error('alpha must be a vector'); +end + +% ensure x is is a vector of the same size as alpha +if any(size(x) ~= size(alpha)) + error('x and alpha must be the same size'); +end + + +%------------------------------------------------- +% Main +%------------------------------------------------- +if any(x < 0) + p = 0; +elseif sum(x) ~= 1 + disp(['dirichletpdf warning: sum(x)~=1, but this may be ' ... + 'due to numerical issues']); + p = 0; +else + z = gammaln(sum(alpha)) - sum(gammaln(alpha)); + z = exp(z); + + p = z * prod(x.^(alpha-1)); +end \ No newline at end of file diff --git a/sourcecodes/bnt-master/KPMstats/dirichletrnd.m b/sourcecodes/bnt-master/KPMstats/dirichletrnd.m new file mode 100644 index 00000000..f11074db --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/dirichletrnd.m @@ -0,0 +1,32 @@ +function x = dirichletrnd(alpha) +%DIRICHLETRND Random vector from a dirichlet distribution. +% x = dirichletrnd(alpha) returns a vector randomly selected +% from the Dirichlet distribution with parameter vector alpha. +% +% The algorithm used is the following: +% For each alpha(i), generate a value s(i) with distribution +% Gamma(alpha(i),1). Now x(i) = s(i) / sum_j s(j). +% +% The above algorithm was recounted to me by Radford Neal, but +% a reference would be appreciated... +% Do the gamma parameters always have to be 1? +% +% Author: David Ross +% $Id: dirichletrnd.m,v 1.1.1.1 2005/05/22 23:32:12 yozhik Exp $ + +%------------------------------------------------- +% Check the input +%------------------------------------------------- +error(nargchk(1,1,nargin)); + +if min(size(alpha)) ~= 1 | length(alpha) < 2 + error('alpha must be a vector of length at least 2'); +end + + +%------------------------------------------------- +% Main +%------------------------------------------------- +gamma_vals = gamrnd(alpha, ones(size(alpha)), size(alpha)); +denom = sum(gamma_vals); +x = gamma_vals / denom; \ No newline at end of file diff --git a/sourcecodes/bnt-master/KPMstats/distchck.m b/sourcecodes/bnt-master/KPMstats/distchck.m new file mode 100644 index 00000000..399f2b3f --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/distchck.m @@ -0,0 +1,173 @@ +function [errorcode,out1,out2,out3,out4] = distchck(nparms,arg1,arg2,arg3,arg4) +%DISTCHCK Checks the argument list for the probability functions. + +% B.A. Jones 1-22-93 +% Copyright (c) 1993-98 by The MathWorks, Inc. +% $Revision: 1.1.1.1 $ $Date: 2005/04/26 02:29:18 $ + +errorcode = 0; + +if nparms == 1 + out1 = arg1; + return; +end + +if nparms == 2 + [r1 c1] = size(arg1); + [r2 c2] = size(arg2); + scalararg1 = (prod(size(arg1)) == 1); + scalararg2 = (prod(size(arg2)) == 1); + if ~scalararg1 & ~scalararg2 + if r1 ~= r2 | c1 ~= c2 + errorcode = 1; + return; + end + end + if scalararg1 + out1 = arg1(ones(r2,1),ones(c2,1)); + else + out1 = arg1; + end + if scalararg2 + out2 = arg2(ones(r1,1),ones(c1,1)); + else + out2 = arg2; + end +end + +if nparms == 3 + [r1 c1] = size(arg1); + [r2 c2] = size(arg2); + [r3 c3] = size(arg3); + scalararg1 = (prod(size(arg1)) == 1); + scalararg2 = (prod(size(arg2)) == 1); + scalararg3 = (prod(size(arg3)) == 1); + + if ~scalararg1 & ~scalararg2 + if r1 ~= r2 | c1 ~= c2 + errorcode = 1; + return; + end + end + + if ~scalararg1 & ~scalararg3 + if r1 ~= r3 | c1 ~= c3 + errorcode = 1; + return; + end + end + + if ~scalararg3 & ~scalararg2 + if r3 ~= r2 | c3 ~= c2 + errorcode = 1; + return; + end + end + + if ~scalararg1 + out1 = arg1; + end + if ~scalararg2 + out2 = arg2; + end + if ~scalararg3 + out3 = arg3; + end + rows = max([r1 r2 r3]); + columns = max([c1 c2 c3]); + + if scalararg1 + out1 = arg1(ones(rows,1),ones(columns,1)); + end + if scalararg2 + out2 = arg2(ones(rows,1),ones(columns,1)); + end + if scalararg3 + out3 = arg3(ones(rows,1),ones(columns,1)); + end + out4 =[]; + +end + +if nparms == 4 + [r1 c1] = size(arg1); + [r2 c2] = size(arg2); + [r3 c3] = size(arg3); + [r4 c4] = size(arg4); + scalararg1 = (prod(size(arg1)) == 1); + scalararg2 = (prod(size(arg2)) == 1); + scalararg3 = (prod(size(arg3)) == 1); + scalararg4 = (prod(size(arg4)) == 1); + + if ~scalararg1 & ~scalararg2 + if r1 ~= r2 | c1 ~= c2 + errorcode = 1; + return; + end + end + + if ~scalararg1 & ~scalararg3 + if r1 ~= r3 | c1 ~= c3 + errorcode = 1; + return; + end + end + + if ~scalararg1 & ~scalararg4 + if r1 ~= r4 | c1 ~= c4 + errorcode = 1; + return; + end + end + + if ~scalararg3 & ~scalararg2 + if r3 ~= r2 | c3 ~= c2 + errorcode = 1; + return; + end + end + + if ~scalararg4 & ~scalararg2 + if r4 ~= r2 | c4 ~= c2 + errorcode = 1; + return; + end + end + + if ~scalararg3 & ~scalararg4 + if r3 ~= r4 | c3 ~= c4 + errorcode = 1; + return; + end + end + + + if ~scalararg1 + out1 = arg1; + end + if ~scalararg2 + out2 = arg2; + end + if ~scalararg3 + out3 = arg3; + end + if ~scalararg4 + out4 = arg4; + end + + rows = max([r1 r2 r3 r4]); + columns = max([c1 c2 c3 c4]); + if scalararg1 + out1 = arg1(ones(rows,1),ones(columns,1)); + end + if scalararg2 + out2 = arg2(ones(rows,1),ones(columns,1)); + end + if scalararg3 + out3 = arg3(ones(rows,1),ones(columns,1)); + end + if scalararg4 + out4 = arg4(ones(rows,1),ones(columns,1)); + end +end + diff --git a/sourcecodes/bnt-master/KPMstats/eigdec.m b/sourcecodes/bnt-master/KPMstats/eigdec.m new file mode 100644 index 00000000..321cee4d --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/eigdec.m @@ -0,0 +1,59 @@ +function [evals, evec] = eigdec(x, N) +%EIGDEC Sorted eigendecomposition +% +% Description +% EVALS = EIGDEC(X, N computes the largest N eigenvalues of the +% matrix X in descending order. [EVALS, EVEC] = EIGDEC(X, N) also +% computes the corresponding eigenvectors. +% +% See also +% PCA, PPCA +% + +% Copyright (c) Ian T Nabney (1996-2001) + +if nargout == 1 + evals_only = logical(1); +else + evals_only = logical(0); +end + +if N ~= round(N) | N < 1 | N > size(x, 2) + error('Number of PCs must be integer, >0, < dim'); +end + +% Find the eigenvalues of the data covariance matrix +if evals_only + % Use eig function as always more efficient than eigs here + temp_evals = eig(x); +else + % Use eig function unless fraction of eigenvalues required is tiny + if (N/size(x, 2)) > 0.04 + fprintf('netlab pca: using eig\n'); + [temp_evec, temp_evals] = eig(x); + else + options.disp = 0; + fprintf('netlab pca: using eigs\n'); + [temp_evec, temp_evals] = eigs(x, N, 'LM', options); + end + temp_evals = diag(temp_evals); +end + +% Eigenvalues nearly always returned in descending order, but just +% to make sure..... +[evals perm] = sort(-temp_evals); +evals = -evals(1:N); +%evec=temp_evec(:,1:N); +if ~evals_only + if evals == temp_evals(1:N) + % Originals were in order + evec = temp_evec(:, 1:N); + return + else + fprintf('netlab pca: sorting evec\n'); + % Need to reorder the eigenvectors + for i=1:N + evec(:,i) = temp_evec(:,perm(i)); + end + end +end diff --git a/sourcecodes/bnt-master/KPMstats/est_transmat.m b/sourcecodes/bnt-master/KPMstats/est_transmat.m new file mode 100644 index 00000000..7654d9ab --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/est_transmat.m @@ -0,0 +1,15 @@ +function [A,C] = est_transmat(seq) +% ESTIMATE_TRANSMAT Max likelihood of a Markov chain transition matrix +% [A,C] = estimate_transmat(seq) +% +% seq is a vector of positive integers +% +% e.g., seq = [1 2 1 2 3], C(1,2)=2, C(2,1)=1, C(2,3)=1, so +% A(1,:)=[0 1 0], A(2,:) = [0.5 0 0.5], +% all other entries are 0 + +% Use a trick with sparse matrices to count the number of each transition. +% From http://www.mathworks.com/company/newsletter/may03/dna.shtml + +C = full(sparse(seq(1:end-1), seq(2:end),1)); +A = mk_stochastic(C); diff --git a/sourcecodes/bnt-master/KPMstats/fit_paritioned_model_testfn.m b/sourcecodes/bnt-master/KPMstats/fit_paritioned_model_testfn.m new file mode 100644 index 00000000..756cdcf9 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/fit_paritioned_model_testfn.m @@ -0,0 +1,5 @@ +function model = foo(inputs, outputs, varargin) + +model.inputs = inputs; +model.outputs = outputs; + diff --git a/sourcecodes/bnt-master/KPMstats/fit_partitioned_model.m b/sourcecodes/bnt-master/KPMstats/fit_partitioned_model.m new file mode 100644 index 00000000..fb9d178d --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/fit_partitioned_model.m @@ -0,0 +1,59 @@ +function [model, partition_size] = fit_partitioned_model(... + inputs, outputs, selectors, sel_sizes, min_size, partition_names, fn_name, varargin) +%function [models, partition_sizes] = fit_partitioned_model(... +% inputs, outputs, selectors, sel_sizes, min_size, partition_names, fn_name, varargin) +% +% Fit models to different subsets (columns) of the input/output data, +% as chosen by the selectors matrix. If there is only output data, set input=[]. +% If there is less than min_size data in partition i, +% we set model{i} = [] +% +% Example: +% selectors = [1 2 1 1 1 +% 1 2 2 1 2] +% sel_sizes = [2 2] so there are 4 models: (1,1), (2,1), (1,2), (2,2) +% We fit model{1} to data from columns 1,4 +% We fit model{2} to no data +% We fit model{3} to data from column 3,5 +% We fit model{4} to data from column 2 (assuming min_size <= 1) +% +% For each partition, we call the specified function with the specified arguments +% as follows: +% model{i} = fn(input(:,cols{i}), output(:,cols{i}), args) +% (We omit input if []) +% partition_size(i) is the amount of data in the i'th partition. +% +% Example use: row 1 of selectors is whether an object is present/absent +% and row 2 is the location. +% +% Demo: +% inputs = 1:5; outputs = 6:10; selectors = as above +% fn = 'fit_partitioned_model_testfn'; +% [model, partition_size] = fit_partitioned_model(inputs, outputs, selectors, [2 2], fn) +% should produce +% model{1}.input = [1 4], model{1}.output = [6 9] +% model{2} = [] +% model{3}.input = [3 5], model{3}.output = [8 10], +% model{4}.input = [2], model{3}.output = [7], +% partition_size = [2 0 2 1] + + +sel_ndx = subv2ind(sel_sizes, selectors'); +Nmodels = prod(sel_sizes); +model = cell(1, Nmodels); +partition_size = zeros(1, Nmodels); +for m=1:Nmodels + ndx = find(sel_ndx==m); + partition_size(m) = length(ndx); + if ~isempty(partition_names) % & (partition_size(m) < min_size) + fprintf('partition %s has size %d, min size = %d\n', ... + partition_names{m}, partition_size(m), min_size); + end + if partition_size(m) >= min_size + if isempty(inputs) + model{m} = feval(fn_name, outputs(:, ndx), varargin{:}); + else + model{m} = feval(fn_name, inputs(:,ndx), outputs(:, ndx), varargin{:}); + end + end +end diff --git a/sourcecodes/bnt-master/KPMstats/gamma_sample.m b/sourcecodes/bnt-master/KPMstats/gamma_sample.m new file mode 100644 index 00000000..e622df40 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/gamma_sample.m @@ -0,0 +1,126 @@ +function r = gamrnd(a,b,m,n); +%GAMRND Random matrices from gamma distribution. +% R = GAMRND(A,B) returns a matrix of random numbers chosen +% from the gamma distribution with parameters A and B. +% The size of R is the common size of A and B if both are matrices. +% If either parameter is a scalar, the size of R is the size of the other +% parameter. Alternatively, R = GAMRND(A,B,M,N) returns an M by N matrix. +% +% Some references refer to the gamma distribution +% with a single parameter. This corresponds to GAMRND +% with B = 1. (See Devroye, pages 401-402.) + +% GAMRND uses a rejection or an inversion method depending on the +% value of A. + +% References: +% [1] L. Devroye, "Non-Uniform Random Variate Generation", +% Springer-Verlag, 1986 + +% B.A. Jones 2-1-93 +% Copyright (c) 1993-98 by The MathWorks, Inc. +% $Revision: 1.1.1.1 $ $Date: 2005/04/26 02:29:18 $ + +if nargin < 2, + error('Requires at least two input arguments.'); +end + + +if nargin == 2 + [errorcode rows columns] = rndcheck(2,2,a,b); +end + +if nargin == 3 + [errorcode rows columns] = rndcheck(3,2,a,b,m); +end + +if nargin == 4 + [errorcode rows columns] = rndcheck(4,2,a,b,m,n); +end + +if errorcode > 0 + error('Size information is inconsistent.'); +end + +% Initialize r to zero. +lth = rows*columns; +r = zeros(lth,1); +a = a(:); b = b(:); + +scalara = (length(a) == 1); +if scalara + a = a*ones(lth,1); +end + +scalarb = (length(b) == 1); +if scalarb + b = b*ones(lth,1); +end + +% If a == 1, then gamma is exponential. (Devroye, page 405). +k = find(a == 1); +if any(k) + r(k) = -b(k) .* log(rand(size(k))); +end + + +k = find(a < 1 & a > 0); +% (Devroye, page 418 Johnk's generator) +if any(k) + c = zeros(lth,1); + d = zeros(lth,1); + c(k) = 1 ./ a(k); + d(k) = 1 ./ (1 - a(k)); + accept = k; + while ~isempty(accept) + u = rand(size(accept)); + v = rand(size(accept)); + x = u .^ c(accept); + y = v .^ d(accept); + k1 = find((x + y) <= 1); + if ~isempty(k1) + e = -log(rand(size(k1))); + r(accept(k1)) = e .* x(k1) ./ (x(k1) + y(k1)); + accept(k1) = []; + end + end + r(k) = r(k) .* b(k); +end + +% Use a rejection method for a > 1. +k = find(a > 1); +% (Devroye, page 410 Best's algorithm) +bb = zeros(size(a)); +c = bb; +if any(k) + bb(k) = a(k) - 1; + c(k) = 3 * a(k) - 3/4; + accept = k; + count = 1; + while ~isempty(accept) + m = length(accept); + u = rand(m,1); + v = rand(m,1); + w = u .* (1 - u); + y = sqrt(c(accept) ./ w) .* (u - 0.5); + x = bb(accept) + y; + k1 = find(x >= 0); + if ~isempty(k1) + z = 64 * (w .^ 3) .* (v .^ 2); + k2 = (z(k1) <= (1 - 2 * (y(k1) .^2) ./ x(k1))); + k3 = k1(find(k2)); + r(accept(k3)) = x(k3); + k4 = k1(find(~k2)); + k5 = k4(find(log(z(k4)) <= (2*(bb(accept(k4)).*log(x(k4)./bb(accept(k4)))-y(k4))))); + r(accept(k5)) = x(k5); + omit = [k3; k5]; + accept(omit) = []; + end + end + r(k) = r(k) .* b(k); +end + +% Return NaN if a or b is not positive. +r(b <= 0 | a <= 0) = NaN; + +r = reshape(r,rows,columns); diff --git a/sourcecodes/bnt-master/KPMstats/gaussian_prob.m b/sourcecodes/bnt-master/KPMstats/gaussian_prob.m new file mode 100644 index 00000000..0ac0b765 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/gaussian_prob.m @@ -0,0 +1,28 @@ +function p = gaussian_prob(x, m, C, use_log) +% GAUSSIAN_PROB Evaluate a multivariate Gaussian density. +% p = gaussian_prob(X, m, C) +% p(i) = N(X(:,i), m, C) where C = covariance matrix and each COLUMN of x is a datavector + +% p = gaussian_prob(X, m, C, 1) returns log N(X(:,i), m, C) (to prevents underflow). +% +% If X has size dxN, then p has size Nx1, where N = number of examples + +if nargin < 4, use_log = 0; end + +if length(m)==1 % scalar + x = x(:)'; +end +[d N] = size(x); +%assert(length(m)==d); % slow +m = m(:); +M = m*ones(1,N); % replicate the mean across columns +denom = (2*pi)^(d/2)*sqrt(abs(det(C))); +mahal = sum(((x-M)'*inv(C)).*(x-M)',2); % Chris Bregler's trick +if any(mahal<0) + warning('mahal < 0 => C is not psd') +end +if use_log + p = -0.5*mahal - log(denom); +else + p = exp(-0.5*mahal) / (denom+eps); +end diff --git a/sourcecodes/bnt-master/KPMstats/gaussian_sample.m b/sourcecodes/bnt-master/KPMstats/gaussian_sample.m new file mode 100644 index 00000000..a2f09df8 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/gaussian_sample.m @@ -0,0 +1,25 @@ +function x = gsamp(mu, covar, nsamp) +%GSAMP Sample from a Gaussian distribution. +% +% Description +% +% X = GSAMP(MU, COVAR, NSAMP) generates a sample of size NSAMP from a +% D-dimensional Gaussian distribution. The Gaussian density has mean +% vector MU and covariance matrix COVAR, and the matrix X has NSAMP +% rows in which each row represents a D-dimensional sample vector. +% +% See also +% GAUSS, DEMGAUSS +% + +% Copyright (c) Ian T Nabney (1996-2001) + +d = size(covar, 1); + +mu = reshape(mu, 1, d); % Ensure that mu is a row vector + +[evec, eval] = eig(covar); + +coeffs = randn(nsamp, d)*sqrt(eval); + +x = ones(nsamp, 1)*mu + coeffs*evec'; diff --git a/sourcecodes/bnt-master/KPMstats/histCmpChi2.m b/sourcecodes/bnt-master/KPMstats/histCmpChi2.m new file mode 100644 index 00000000..58f3da98 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/histCmpChi2.m @@ -0,0 +1,15 @@ +function d = histCmpChi2(h1, h2) +% Compare two histograms using chi-squared +% function d = histCmpChi2(h1, h2) +% +% d(i,j) = chi^2(h1(i,:), h2(j,:)) = sum_b (h1(i,b)-h2(j,b)^2 / (h1(i,b) + h2(j,b)) + +[N B] = size(h1); +d = zeros(N,N); +for i=1:N + h1i = repmat(h1(i,:), N, 1); + numer = (h1i - h2).^2; + denom = h1i + h2 + eps; % if denom=0, then numer=0 + d(i,:) = sum(numer ./ denom, 2); +end + diff --git a/sourcecodes/bnt-master/KPMstats/linear_regression.m b/sourcecodes/bnt-master/KPMstats/linear_regression.m new file mode 100644 index 00000000..b33ab774 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/linear_regression.m @@ -0,0 +1,68 @@ +function [muY, SigmaY, weightsY] = linear_regression(X, Y, varargin) +% LINEAR_REGRESSION Fit params for P(Y|X) = N(Y; W X + mu, Sigma) +% +% X(:, t) is the t'th input example +% Y(:, t) is the t'th output example +% +% Kevin Murphy, August 2003 +% +% This is a special case of cwr_em with 1 cluster. +% You can also think of it as a front end to clg_Mstep. + +[cov_typeY, clamp_weights, muY, SigmaY, weightsY,... + cov_priorY, regress, clamp_covY] = process_options(... + varargin, ... + 'cov_typeY', 'full', 'clamp_weights', 0, ... + 'muY', [], 'SigmaY', [], 'weightsY', [], ... + 'cov_priorY', [], 'regress', 1, 'clamp_covY', 0); + +[nx N] = size(X); +[ny N2] = size(Y); +if N ~= N2 + error(sprintf('nsamples X (%d) ~= nsamples Y (%d)', N, N2)); +end + +w = 1/N; +WYbig = Y*w; +WYY = WYbig * Y'; +WY = sum(WYbig, 2); +WYTY = sum(diag(WYbig' * Y)); +if ~regress + % This is just fitting an unconditional Gaussian + weightsY = []; + [muY, SigmaY] = ... + mixgauss_Mstep(1, WY, WYY, WYTY, ... + 'cov_type', cov_typeY, 'cov_prior', cov_priorY); + % There is a much easier way... + assert(approxeq(muY, mean(Y'))) + assert(approxeq(SigmaY, cov(Y') + 0.01*eye(ny))) +else + % This is just linear regression + WXbig = X*w; + WXX = WXbig * X'; + WX = sum(WXbig, 2); + WXTX = sum(diag(WXbig' * X)); + WXY = WXbig * Y'; + [muY, SigmaY, weightsY] = ... + clg_Mstep(1, WY, WYY, WYTY, WX, WXX, WXY, ... + 'cov_type', cov_typeY, 'cov_prior', cov_priorY); +end +if clamp_covY, SigmaY = SigmaY; end +if clamp_weights, weightsY = weightsY; end + +if nx==1 & ny==1 & regress + P = polyfit(X,Y); % Y = P(1) X^1 + P(2) X^0 = ax + b + assert(approxeq(muY, P(2))) + assert(approxeq(weightsY, P(1))) +end + +%%%%%%%% Test +if 0 + c1 = randn(2,100); c2 = randn(2,100); + y = c2(1,:); X = [ones(size(c1,2),1) c1']; + b = regress(y(:), X); % stats toolbox + [m,s,w] = linear_regression(c1, y); + assert(approxeq(b(1),m)) + assert(approxeq(b(2), w(1))) + assert(approxeq(b(3), w(2))) +end diff --git a/sourcecodes/bnt-master/KPMstats/logist2.m b/sourcecodes/bnt-master/KPMstats/logist2.m new file mode 100644 index 00000000..1c1e1619 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/logist2.m @@ -0,0 +1,115 @@ +function [beta,p,lli] = logist2(y,x,w) +% [beta,p,lli] = logist2(y,x) +% +% 2-class logistic regression. +% +% INPUT +% y Nx1 colum vector of 0|1 class assignments +% x NxK matrix of input vectors as rows +% [w] Nx1 vector of sample weights +% +% OUTPUT +% beta Kx1 column vector of model coefficients +% p Nx1 column vector of fitted class 1 posteriors +% lli log likelihood +% +% Class 1 posterior is 1 / (1 + exp(-x*beta)) +% +% David Martin <dmartin@eecs.berkeley.edu> +% April 16, 2002 + +% Copyright (C) 2002 David R. Martin <dmartin@eecs.berkeley.edu> +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. +% + +% This program is distributed in the hope that it will be useful, but +% WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +% General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program; if not, write to the Free Software +% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA +% 02111-1307, USA, or see http://www.gnu.org/copyleft/gpl.html. + +error(nargchk(2,3,nargin)); + +% check inputs +if size(y,2) ~= 1, + error('Input y not a column vector.'); +end +if size(y,1) ~= size(x,1), + error('Input x,y sizes mismatched.'); +end + +% get sizes +[N,k] = size(x); + +% if sample weights weren't specified, set them to 1 +if nargin < 3, + w = 1; +end + +% normalize sample weights so max is 1 +w = w / max(w); + +% initial guess for beta: all zeros +beta = zeros(k,1); + +% Newton-Raphson via IRLS, +% taken from Hastie/Tibshirani/Friedman Section 4.4. +iter = 0; +lli = 0; +while 1==1, + iter = iter + 1; + + % fitted probabilities + p = 1 ./ (1 + exp(-x*beta)); + + % log likelihood + lli_prev = lli; + lli = sum( w .* (y.*log(p+eps) + (1-y).*log(1-p+eps)) ); + + % least-squares weights + wt = w .* p .* (1-p); + + % derivatives of likelihood w.r.t. beta + deriv = x'*(w.*(y-p)); + + % Hessian of likelihood w.r.t. beta + % hessian = x'Wx, where W=diag(w) + % Do it this way to be memory efficient and fast. + hess = zeros(k,k); + for i = 1:k, + wxi = wt .* x(:,i); + for j = i:k, + hij = wxi' * x(:,j); + hess(i,j) = -hij; + hess(j,i) = -hij; + end + end + + % make sure Hessian is well conditioned + if (rcond(hess) < eps), + error(['Stopped at iteration ' num2str(iter) ... + ' because Hessian is poorly conditioned.']); + break; + end; + + % Newton-Raphson update step + step = hess\deriv; + beta = beta - step; + + % termination criterion based on derivatives + tol = 1e-6; + if abs(deriv'*step/k) < tol, break; end; + + % termination criterion based on log likelihood +% tol = 1e-4; +% if abs((lli-lli_prev)/(lli+lli_prev)) < 0.5*tol, break; end; +end; + diff --git a/sourcecodes/bnt-master/KPMstats/logist2Apply.m b/sourcecodes/bnt-master/KPMstats/logist2Apply.m new file mode 100644 index 00000000..57be8c35 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/logist2Apply.m @@ -0,0 +1,13 @@ +function p = logist2Apply(beta, x) +% LOGIST2APPLY 2 class logistic regression: compute posterior prob of class 1 +% function p = logist2Apply(beta, x) +% +% x(:,i) - each COLUMN is a test case; we append 1s automatically, if appropriate + +[D Ncases] = size(x); +if length(beta)==D+1 + F = [x; ones(1,Ncases)]; +else + F = x; +end +p = 1./(1+exp(-beta(:)'*F)); diff --git a/sourcecodes/bnt-master/KPMstats/logist2ApplyRegularized.m b/sourcecodes/bnt-master/KPMstats/logist2ApplyRegularized.m new file mode 100644 index 00000000..39095a3d --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/logist2ApplyRegularized.m @@ -0,0 +1,3 @@ +function prob = logist2ApplyRegularized(net, features) + +prob = glmfwd(net, features')'; diff --git a/sourcecodes/bnt-master/KPMstats/logist2Fit.m b/sourcecodes/bnt-master/KPMstats/logist2Fit.m new file mode 100644 index 00000000..7bc466b7 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/logist2Fit.m @@ -0,0 +1,22 @@ +function [beta, p] = logist2Fit(y, x, addOne, w) +% LOGIST2FIT 2 class logsitic classification +% function beta = logist2Fit(y,x, addOne) +% +% y(i) = 0/1 +% x(:,i) = i'th input - we optionally append 1s to last dimension +% w(i) = optional weight +% +% beta(j)- regression coefficient + +if nargin < 3, addOne = 1; end +if nargin < 4, w = 1; end + +Ncases = size(x,2); +if Ncases ~= length(y) + error(sprintf('size of data = %dx%d, size of labels=%d', size(x,1), size(x,2), length(y))) +end +if addOne + x = [x; ones(1,Ncases)]; +end +[beta, p] = logist2(y(:), x', w(:)); +beta = beta(:); diff --git a/sourcecodes/bnt-master/KPMstats/logist2FitRegularized.m b/sourcecodes/bnt-master/KPMstats/logist2FitRegularized.m new file mode 100644 index 00000000..d521e8f7 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/logist2FitRegularized.m @@ -0,0 +1,13 @@ +function [net, niter] = logist2FitRegularized(labels, features, maxIter) + +if nargin < 3, maxIter = 100; end + +[D N] = size(features); +weightPrior = 0.5; +net = glm(D, 1, 'logistic', weightPrior); +options = foptions; +options(14) = maxIter; +[net, options] = glmtrain(net, options, features', labels(:)); +niter = options(14); +%w = logist2Fit(labelsPatches(jValidPatches), features(:, jValidPatches)); + diff --git a/sourcecodes/bnt-master/KPMstats/logistK.m b/sourcecodes/bnt-master/KPMstats/logistK.m new file mode 100644 index 00000000..8309328a --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/logistK.m @@ -0,0 +1,287 @@ +function [beta,post,lli] = logistK(x,y,w,beta) +% [beta,post,lli] = logistK(x,y,beta,w) +% +% k-class logistic regression with optional sample weights +% +% k = number of classes +% n = number of samples +% d = dimensionality of samples +% +% INPUT +% x dxn matrix of n input column vectors +% y kxn vector of class assignments +% [w] 1xn vector of sample weights +% [beta] dxk matrix of model coefficients +% +% OUTPUT +% beta dxk matrix of fitted model coefficients +% (beta(:,k) are fixed at 0) +% post kxn matrix of fitted class posteriors +% lli log likelihood +% +% Let p(i,j) = exp(beta(:,j)'*x(:,i)), +% Class j posterior for observation i is: +% post(j,i) = p(i,j) / (p(i,1) + ... p(i,k)) +% +% See also logistK_eval. +% +% David Martin <dmartin@eecs.berkeley.edu> +% May 3, 2002 + +% Copyright (C) 2002 David R. Martin <dmartin@eecs.berkeley.edu> +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, but +% WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +% General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program; if not, write to the Free Software +% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA +% 02111-1307, USA, or see http://www.gnu.org/copyleft/gpl.html. + +% TODO - this code would be faster if x were transposed + +error(nargchk(2,4,nargin)); + +debug = 0; +if debug>0, + h=figure(1); + set(h,'DoubleBuffer','on'); +end + +% get sizes +[d,nx] = size(x); +[k,ny] = size(y); + +% check sizes +if k < 2, + error('Input y must encode at least 2 classes.'); +end +if nx ~= ny, + error('Inputs x,y not the same length.'); +end + +n = nx; + +% make sure class assignments have unit L1-norm +sumy = sum(y,1); +if abs(1-sumy) > eps, + sumy = sum(y,1); + for i = 1:k, y(i,:) = y(i,:) ./ sumy; end +end +clear sumy; + +% if sample weights weren't specified, set them to 1 +if nargin < 3, + w = ones(1,n); +end + +% normalize sample weights so max is 1 +w = w / max(w); + +% if starting beta wasn't specified, initialize randomly +if nargin < 4, + beta = 1e-3*rand(d,k); + beta(:,k) = 0; % fix beta for class k at zero +else + if sum(beta(:,k)) ~= 0, + error('beta(:,k) ~= 0'); + end +end + +stepsize = 1; +minstepsize = 1e-2; + +post = computePost(beta,x); +lli = computeLogLik(post,y,w); + +for iter = 1:100, + %disp(sprintf(' logist iter=%d lli=%g',iter,lli)); + vis(x,y,beta,lli,d,k,iter,debug); + + % gradient and hessian + [g,h] = derivs(post,x,y,w); + + % make sure Hessian is well conditioned + if rcond(h) < eps, + % condition with Levenberg-Marquardt method + for i = -16:16, + h2 = h .* ((1 + 10^i)*eye(size(h)) + (1-eye(size(h)))); + if rcond(h2) > eps, break, end + end + if rcond(h2) < eps, + warning(['Stopped at iteration ' num2str(iter) ... + ' because Hessian can''t be conditioned']); + break + end + h = h2; + end + + % save lli before update + lli_prev = lli; + + % Newton-Raphson with step-size halving + while stepsize >= minstepsize, + % Newton-Raphson update step + step = stepsize * (h \ g); + beta2 = beta; + beta2(:,1:k-1) = beta2(:,1:k-1) - reshape(step,d,k-1); + + % get the new log likelihood + post2 = computePost(beta2,x); + lli2 = computeLogLik(post2,y,w); + + % if the log likelihood increased, then stop + if lli2 > lli, + post = post2; lli = lli2; beta = beta2; + break + end + + % otherwise, reduce step size by half + stepsize = 0.5 * stepsize; + end + + % stop if the average log likelihood has gotten small enough + if 1-exp(lli/n) < 1e-2, break, end + + % stop if the log likelihood changed by a small enough fraction + dlli = (lli_prev-lli) / lli; + if abs(dlli) < 1e-3, break, end + + % stop if the step size has gotten too small + if stepsize < minstepsize, brea, end + + % stop if the log likelihood has decreased; this shouldn't happen + if lli < lli_prev, + warning(['Stopped at iteration ' num2str(iter) ... + ' because the log likelihood decreased from ' ... + num2str(lli_prev) ' to ' num2str(lli) '.' ... + ' This may be a bug.']); + break + end +end + +if debug>0, + vis(x,y,beta,lli,d,k,iter,2); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% class posteriors +function post = computePost(beta,x) + [d,n] = size(x); + [d,k] = size(beta); + post = zeros(k,n); + bx = zeros(k,n); + for j = 1:k, + bx(j,:) = beta(:,j)'*x; + end + for j = 1:k, + post(j,:) = 1 ./ sum(exp(bx - repmat(bx(j,:),k,1)),1); + end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% log likelihood +function lli = computeLogLik(post,y,w) + [k,n] = size(post); + lli = 0; + for j = 1:k, + lli = lli + sum(w.*y(j,:).*log(post(j,:)+eps)); + end + if isnan(lli), + error('lli is nan'); + end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% gradient and hessian +%% These are computed in what seems a verbose manner, but it is +%% done this way to use minimal memory. x should be transposed +%% to make it faster. +function [g,h] = derivs(post,x,y,w) + + [k,n] = size(post); + [d,n] = size(x); + + % first derivative of likelihood w.r.t. beta + g = zeros(d,k-1); + for j = 1:k-1, + wyp = w .* (y(j,:) - post(j,:)); + for ii = 1:d, + g(ii,j) = x(ii,:) * wyp'; + end + end + g = reshape(g,d*(k-1),1); + + % hessian of likelihood w.r.t. beta + h = zeros(d*(k-1),d*(k-1)); + for i = 1:k-1, % diagonal + wt = w .* post(i,:) .* (1 - post(i,:)); + hii = zeros(d,d); + for a = 1:d, + wxa = wt .* x(a,:); + for b = a:d, + hii_ab = wxa * x(b,:)'; + hii(a,b) = hii_ab; + hii(b,a) = hii_ab; + end + end + h( (i-1)*d+1 : i*d , (i-1)*d+1 : i*d ) = -hii; + end + for i = 1:k-1, % off-diagonal + for j = i+1:k-1, + wt = w .* post(j,:) .* post(i,:); + hij = zeros(d,d); + for a = 1:d, + wxa = wt .* x(a,:); + for b = a:d, + hij_ab = wxa * x(b,:)'; + hij(a,b) = hij_ab; + hij(b,a) = hij_ab; + end + end + h( (i-1)*d+1 : i*d , (j-1)*d+1 : j*d ) = hij; + h( (j-1)*d+1 : j*d , (i-1)*d+1 : i*d ) = hij; + end + end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% debug/visualization +function vis (x,y,beta,lli,d,k,iter,debug) + + if debug<=0, return, end + + disp(['iter=' num2str(iter) ' lli=' num2str(lli)]); + if debug<=1, return, end + + if d~=3 | k>10, return, end + + figure(1); + res = 100; + r = abs(max(max(x))); + dom = linspace(-r,r,res); + [px,py] = meshgrid(dom,dom); + xx = px(:); yy = py(:); + points = [xx' ; yy' ; ones(1,res*res)]; + func = zeros(k,res*res); + for j = 1:k, + func(j,:) = exp(beta(:,j)'*points); + end + [mval,ind] = max(func,[],1); + hold off; + im = reshape(ind,res,res); + imagesc(xx,yy,im); + hold on; + syms = {'w.' 'wx' 'w+' 'wo' 'w*' 'ws' 'wd' 'wv' 'w^' 'w<'}; + for j = 1:k, + [mval,ind] = max(y,[],1); + ind = find(ind==j); + plot(x(1,ind),x(2,ind),syms{j}); + end + pause(0.1); + +% eof diff --git a/sourcecodes/bnt-master/KPMstats/logistK_eval.m b/sourcecodes/bnt-master/KPMstats/logistK_eval.m new file mode 100644 index 00000000..880a5f4b --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/logistK_eval.m @@ -0,0 +1,83 @@ +function [post,lik,lli] = logistK_eval(beta,x,y) +% [post,lik,lli] = logistK_eval(beta,x,y) +% +% Evaluate logistic regression model. +% +% INPUT +% beta dxk model coefficients (as returned by logistK) +% x dxn matrix of n input column vectors +% [y] kxn vector of class assignments +% +% OUTPUT +% post kxn fitted class posteriors +% lik 1xn vector of sample likelihoods +% lli log likelihood +% +% Let p(i,j) = exp(beta(:,j)'*x(:,i)), +% Class j posterior for observation i is: +% post(j,i) = p(i,j) / (p(i,1) + ... p(i,k)) +% The likelihood of observation i given soft class assignments +% y(:,i) is: +% lik(i) = prod(post(:,i).^y(:,i)) +% The log-likelihood of the model given the labeled samples is: +% lli = sum(log(lik)) +% +% See also logistK. +% +% David Martin <dmartin@eecs.berkeley.edu> +% May 7, 2002 + +% Copyright (C) 2002 David R. Martin <dmartin@eecs.berkeley.edu> +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, but +% WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +% General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program; if not, write to the Free Software +% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA +% 02111-1307, USA, or see http://www.gnu.org/copyleft/gpl.html. + +error(nargchk(2,3,nargin)); + +% check sizes +if size(beta,1) ~= size(x,1), + error('Inputs beta,x not the same height.'); +end +if nargin > 3 & size(y,2) ~= size(x,2), + error('Inputs x,y not the same length.'); +end + +% get sizes +[d,k] = size(beta); +[d,n] = size(x); + +% class posteriors +post = zeros(k,n); +bx = zeros(k,n); +for j = 1:k, + bx(j,:) = beta(:,j)'*x; +end +for j = 1:k, + post(j,:) = 1 ./ sum(exp(bx - repmat(bx(j,:),k,1)),1); +end +clear bx; + +% likelihood of each sample +if nargout > 1, + y = y ./ repmat(sum(y,1),k,1); % L1-normalize class assignments + lik = prod(post.^y,1); +end + +% total log likelihood +if nargout > 2, + lli = sum(log(lik+eps)); +end; + +% eof diff --git a/sourcecodes/bnt-master/KPMstats/marginalize_gaussian.m b/sourcecodes/bnt-master/KPMstats/marginalize_gaussian.m new file mode 100644 index 00000000..7f347f48 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/marginalize_gaussian.m @@ -0,0 +1,7 @@ +function [muX, SXX] = marginalize_gaussian(mu, Sigma, X, Y, ns) +% MARGINALIZE_GAUSSIAN Compute Pr(X) from Pr(X,Y) where X and Y are jointly Gaussian. +% [muX, SXX] = marginalize_gaussian(mu, Sigma, X, Y, ns) + +[muX, muY, SXX, SXY, SYX, SYY] = partition_matrix_vec(mu, Sigma, X, Y, ns); + + diff --git a/sourcecodes/bnt-master/KPMstats/matrix_T_pdf.m b/sourcecodes/bnt-master/KPMstats/matrix_T_pdf.m new file mode 100644 index 00000000..b6fb92d4 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/matrix_T_pdf.m @@ -0,0 +1,12 @@ +function p = matrix_T_pdf(A, M, V, K, n) +% MATRIX_T_PDF Evaluate the density of a matrix under a Matrix-T distribution +% p = matrix_T_pdf(A, M, V, K, n) + +% See "Bayesian Linear Regression", T. Minka, MIT Tech Report, 2001 + +[d m] = size(K); +is = 1:d; +c1 = prod(gamma((n+1-is)/2)) / prod(gamma((n-m+1-is)/2)); +c2 = det(K)^(d/2) / det(pi*V)^(m/2); %% pi or 2pi? +p = c1 * c2 * det((A-M)'*inv(V)*(A-M)*K + eye(m))^(-n/2); + diff --git a/sourcecodes/bnt-master/KPMstats/matrix_normal_pdf.m b/sourcecodes/bnt-master/KPMstats/matrix_normal_pdf.m new file mode 100644 index 00000000..5e8e9a78 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/matrix_normal_pdf.m @@ -0,0 +1,9 @@ +function p = matrix_normal_pdf(A, M, V, K) +% MATRIX_NORMAL_PDF Evaluate the density of a matrix under a Matrix-Normal distribution +% p = matrix_normal_pdf(A, M, V, K) + +% See "Bayesian Linear Regression", T. Minka, MIT Tech Report, 2001 + +[d m] = size(K); +c = det(K)^(d/2) / det(2*pi*V)^(m/2); +p = c * exp(-0.5*tr((A-M)'*inv(V)*(A-M)*K)); diff --git a/sourcecodes/bnt-master/KPMstats/mc_stat_distrib.m b/sourcecodes/bnt-master/KPMstats/mc_stat_distrib.m new file mode 100644 index 00000000..a5806092 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/mc_stat_distrib.m @@ -0,0 +1,26 @@ +function pi = mc_stat_distrib(P) +% MC_STAT_DISTRIB Compute stationary distribution of a Markov chain +% function pi = mc_stat_distrib(P) +% +% Each row of P should sum to one; pi is a column vector + +% Kevin Murphy, 16 Feb 2003 + +% The stationary distribution pi satisfies pi P = pi +% subject to sum_i pi(i) = 1, 0 <= pi(i) <= 1 +% Hence +% (P' 0n (pi = (pi +% 1n 0) 1) 1) +% or P2 pi2 = pi2. +% Naively we can solve this using (P2 - I(n+1)) pi2 = 0(n+1) +% or P3 pi2 = 0(n+1), i.e., pi2 = P3 \ zeros(n+1,1) +% but this is singular (because of the sum-to-one constraint). +% Hence we replace the last row of P' with 1s instead of appending ones to create P2, +% and similarly for pi. + +n = length(P); +P4 = P'-eye(n); +P4(end,:) = 1; +pi = P4 \ [zeros(n-1,1);1]; + + diff --git a/sourcecodes/bnt-master/KPMstats/mixgauss_Mstep.m b/sourcecodes/bnt-master/KPMstats/mixgauss_Mstep.m new file mode 100644 index 00000000..7908a704 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/mixgauss_Mstep.m @@ -0,0 +1,106 @@ +function [mu, Sigma] = mixgauss_Mstep(w, Y, YY, YTY, varargin) +% MSTEP_COND_GAUSS Compute MLEs for mixture of Gaussians given expected sufficient statistics +% function [mu, Sigma] = Mstep_cond_gauss(w, Y, YY, YTY, varargin) +% +% We assume P(Y|Q=i) = N(Y; mu_i, Sigma_i) +% and w(i,t) = p(Q(t)=i|y(t)) = posterior responsibility +% See www.ai.mit.edu/~murphyk/Papers/learncg.pdf. +% +% INPUTS: +% w(i) = sum_t w(i,t) = responsibilities for each mixture component +% If there is only one mixture component (i.e., Q does not exist), +% then w(i) = N = nsamples, and +% all references to i can be replaced by 1. +% YY(:,:,i) = sum_t w(i,t) y(:,t) y(:,t)' = weighted outer product +% Y(:,i) = sum_t w(i,t) y(:,t) = weighted observations +% YTY(i) = sum_t w(i,t) y(:,t)' y(:,t) = weighted inner product +% You only need to pass in YTY if Sigma is to be estimated as spherical. +% +% Optional parameters may be passed as 'param_name', param_value pairs. +% Parameter names are shown below; default values in [] - if none, argument is mandatory. +% +% 'cov_type' - 'full', 'diag' or 'spherical' ['full'] +% 'tied_cov' - 1 (Sigma) or 0 (Sigma_i) [0] +% 'clamped_cov' - pass in clamped value, or [] if unclamped [ [] ] +% 'clamped_mean' - pass in clamped value, or [] if unclamped [ [] ] +% 'cov_prior' - Lambda_i, added to YY(:,:,i) [0.01*eye(d,d,Q)] +% +% If covariance is tied, Sigma has size d*d. +% But diagonal and spherical covariances are represented in full size. + +[cov_type, tied_cov, clamped_cov, clamped_mean, cov_prior, other] = ... + process_options(varargin,... + 'cov_type', 'full', 'tied_cov', 0, 'clamped_cov', [], 'clamped_mean', [], ... + 'cov_prior', []); + +[Ysz Q] = size(Y); +N = sum(w); +if isempty(cov_prior) + %cov_prior = zeros(Ysz, Ysz, Q); + %for q=1:Q + % cov_prior(:,:,q) = 0.01*cov(Y(:,q)'); + %end + cov_prior = repmat(0.01*eye(Ysz,Ysz), [1 1 Q]); +end +%YY = reshape(YY, [Ysz Ysz Q]) + cov_prior; % regularize the scatter matrix +YY = reshape(YY, [Ysz Ysz Q]); + +% Set any zero weights to one before dividing +% This is valid because w(i)=0 => Y(:,i)=0, etc +w = w + (w==0); + +if ~isempty(clamped_mean) + mu = clamped_mean; +else + % eqn 6 + %mu = Y ./ repmat(w(:)', [Ysz 1]);% Y may have a funny size + mu = zeros(Ysz, Q); + for i=1:Q + mu(:,i) = Y(:,i) / w(i); + end +end + +if ~isempty(clamped_cov) + Sigma = clamped_cov; + return; +end + +if ~tied_cov + Sigma = zeros(Ysz,Ysz,Q); + for i=1:Q + if cov_type(1) == 's' + % eqn 17 + s2 = (1/Ysz)*( (YTY(i)/w(i)) - mu(:,i)'*mu(:,i) ); + Sigma(:,:,i) = s2 * eye(Ysz); + else + % eqn 12 + SS = YY(:,:,i)/w(i) - mu(:,i)*mu(:,i)'; + if cov_type(1)=='d' + SS = diag(diag(SS)); + end + Sigma(:,:,i) = SS; + end + end +else % tied cov + if cov_type(1) == 's' + % eqn 19 + s2 = (1/(N*Ysz))*(sum(YTY,2) + sum(diag(mu'*mu) .* w)); + Sigma = s2*eye(Ysz); + else + SS = zeros(Ysz, Ysz); + % eqn 15 + for i=1:Q % probably could vectorize this... + SS = SS + YY(:,:,i)/N - mu(:,i)*mu(:,i)'; + end + if cov_type(1) == 'd' + Sigma = diag(diag(SS)); + else + Sigma = SS; + end + end +end + +if tied_cov + Sigma = repmat(Sigma, [1 1 Q]); +end +Sigma = Sigma + cov_prior; diff --git a/sourcecodes/bnt-master/KPMstats/mixgauss_classifier_apply.m b/sourcecodes/bnt-master/KPMstats/mixgauss_classifier_apply.m new file mode 100644 index 00000000..73baaeb1 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/mixgauss_classifier_apply.m @@ -0,0 +1,11 @@ +function [classHatTest, probPos] = mixgauss_classifier_apply(mixgauss, testFeatures) + +Bpos = mixgauss_prob(testFeatures, mixgauss.pos.mu, mixgauss.pos.Sigma, mixgauss.pos.prior); +Bneg = mixgauss_prob(testFeatures, mixgauss.neg.mu, mixgauss.neg.Sigma, mixgauss.neg.prior); +prior_pos = mixgauss.priorC(1); +prior_neg = mixgauss.priorC(2); +post = normalize([Bpos * prior_pos; Bneg * prior_neg], 1); +probPos = post(1,:)'; +[junk, classHatTest] = max(post); +classHatTest(find(classHatTest==2))=0; +classHatTest = classHatTest(:); diff --git a/sourcecodes/bnt-master/KPMstats/mixgauss_classifier_train.m b/sourcecodes/bnt-master/KPMstats/mixgauss_classifier_train.m new file mode 100644 index 00000000..fa0861ac --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/mixgauss_classifier_train.m @@ -0,0 +1,33 @@ +function mixgauss = mixgauss_classifier_train(trainFeatures, trainLabels, nc, varargin) +% function mixgauss = mixgauss_classifier_train(trainFeatures, trainLabels, nclusters, varargin) +% trainFeatures(:,i) for i'th example +% trainLabels should be 0,1 +% To evaluate performance on a tets set, use +% mixgauss = mixgauss_classifier_train(trainFeatures, trainLabels, nc, 'testFeatures', tf, 'testLabels', tl) + +[testFeatures, testLabels, max_iter, thresh, cov_type, mu, Sigma, priorC, method, ... + cov_prior, verbose, prune_thresh] = process_options(... + varargin, 'testFeatures', [], 'testLabels', [], ... + 'max_iter', 10, 'thresh', 0.01, 'cov_type', 'diag', ... + 'mu', [], 'Sigma', [], 'priorC', [], 'method', 'kmeans', ... + 'cov_prior', [], 'verbose', 0, 'prune_thresh', 0); + +Nclasses = 2; % max([trainLabels testLabels]) + 1; + +pos = find(trainLabels == 1); +neg = find(trainLabels == 0); + +if verbose, fprintf('fitting pos\n'); end +[mixgauss.pos.mu, mixgauss.pos.Sigma, mixgauss.pos.prior] = ... + mixgauss_em(trainFeatures(:, pos), nc, varargin{:}); + +if verbose, fprintf('fitting neg\n'); end +[mixgauss.neg.mu, mixgauss.neg.Sigma, mixgauss.neg.prior] = ... + mixgauss_em(trainFeatures(:, neg), nc, varargin{:}); + + +if ~isempty(priorC) + mixgauss.priorC = priorC; +else + mixgauss.priorC = normalize([length(pos) length(neg)]); +end diff --git a/sourcecodes/bnt-master/KPMstats/mixgauss_em.m b/sourcecodes/bnt-master/KPMstats/mixgauss_em.m new file mode 100644 index 00000000..a0a8154d --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/mixgauss_em.m @@ -0,0 +1,74 @@ +function [mu, Sigma, prior] = mixgauss_em(Y, nc, varargin) +% MIXGAUSS_EM Fit the parameters of a mixture of Gaussians using EM +% function [mu, Sigma, prior] = mixgauss_em(data, nc, varargin) +% +% data(:, t) is the t'th data point +% nc is the number of clusters + +% Kevin Murphy, 13 May 2003 + +[max_iter, thresh, cov_type, mu, Sigma, method, ... + cov_prior, verbose, prune_thresh] = process_options(... + varargin, 'max_iter', 10, 'thresh', 1e-2, 'cov_type', 'full', ... + 'mu', [], 'Sigma', [], 'method', 'kmeans', ... + 'cov_prior', [], 'verbose', 0, 'prune_thresh', 0); + +[ny T] = size(Y); + +if nc==1 + % No latent variable, so there is a closed-form solution + mu = mean(Y')'; + Sigma = cov(Y'); + if strcmp(cov_type, 'diag') + Sigma = diag(diag(Sigma)); + end + prior = 1; + return; +end + +if isempty(mu) + [mu, Sigma, prior] = mixgauss_init(nc, Y, cov_type, method); +end + +previous_loglik = -inf; +num_iter = 1; +converged = 0; + +%if verbose, fprintf('starting em\n'); end + +while (num_iter <= max_iter) & ~converged + % E step + probY = mixgauss_prob(Y, mu, Sigma, prior); % probY(q,t) + [post, lik] = normalize(probY .* repmat(prior, 1, T), 1); % post(q,t) + loglik = log(sum(lik)); + + % extract expected sufficient statistics + w = sum(post,2); % w(c) = sum_t post(c,t) + WYY = zeros(ny, ny, nc); % WYY(:,:,c) = sum_t post(c,t) Y(:,t) Y(:,t)' + WY = zeros(ny, nc); % WY(:,c) = sum_t post(c,t) Y(:,t) + WYTY = zeros(nc,1); % WYTY(c) = sum_t post(c,t) Y(:,t)' Y(:,t) + for c=1:nc + weights = repmat(post(c,:), ny, 1); % weights(:,t) = post(c,t) + WYbig = Y .* weights; % WYbig(:,t) = post(c,t) * Y(:,t) + WYY(:,:,c) = WYbig * Y'; + WY(:,c) = sum(WYbig, 2); + WYTY(c) = sum(diag(WYbig' * Y)); + end + + % M step + prior = normalize(w); + [mu, Sigma] = mixgauss_Mstep(w, WY, WYY, WYTY, 'cov_type', cov_type, 'cov_prior', cov_prior); + + if verbose, fprintf(1, 'iteration %d, loglik = %f\n', num_iter, loglik); end + num_iter = num_iter + 1; + converged = em_converged(loglik, previous_loglik, thresh); + previous_loglik = loglik; + +end + +if prune_thresh > 0 + ndx = find(prior < prune_thresh); + mu(:,ndx) = []; + Sigma(:,:,ndx) = []; + prior(ndx) = []; +end diff --git a/sourcecodes/bnt-master/KPMstats/mixgauss_init.m b/sourcecodes/bnt-master/KPMstats/mixgauss_init.m new file mode 100644 index 00000000..5a4596e1 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/mixgauss_init.m @@ -0,0 +1,49 @@ +function [mu, Sigma, weights] = mixgauss_init(M, data, cov_type, method) +% MIXGAUSS_INIT Initial parameter estimates for a mixture of Gaussians +% function [mu, Sigma, weights] = mixgauss_init(M, data, cov_type. method) +% +% INPUTS: +% data(:,t) is the t'th example +% M = num. mixture components +% cov_type = 'full', 'diag' or 'spherical' +% method = 'rnd' (choose centers randomly from data) or 'kmeans' (needs netlab) +% +% OUTPUTS: +% mu(:,k) +% Sigma(:,:,k) +% weights(k) + +if nargin < 4, method = 'kmeans'; end + +[d T] = size(data); +data = reshape(data, d, T); % in case it is data(:, t, sequence_num) + +switch method + case 'rnd', + C = cov(data'); + Sigma = repmat(diag(diag(C))*0.5, [1 1 M]); + % Initialize each mean to a random data point + indices = randperm(T); + mu = data(:,indices(1:M)); + weights = normalise(ones(M,1)); + case 'kmeans', + mix = gmm(d, M, cov_type); + options = foptions; + max_iter = 5; + options(1) = -1; % be quiet! + options(14) = max_iter; + mix = gmminit(mix, data', options); + mu = reshape(mix.centres', [d M]); + weights = mix.priors(:); + for m=1:M + switch cov_type + case 'diag', + Sigma(:,:,m) = diag(mix.covars(m,:)); + case 'full', + Sigma(:,:,m) = mix.covars(:,:,m); + case 'spherical', + Sigma(:,:,m) = mix.covars(m) * eye(d); + end + end +end + diff --git a/sourcecodes/bnt-master/KPMstats/mixgauss_prob.m b/sourcecodes/bnt-master/KPMstats/mixgauss_prob.m new file mode 100644 index 00000000..3bd5afcc --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/mixgauss_prob.m @@ -0,0 +1,133 @@ +function [B, B2] = mixgauss_prob(data, mu, Sigma, mixmat, unit_norm) +% EVAL_PDF_COND_MOG Evaluate the pdf of a conditional mixture of Gaussians +% function [B, B2] = eval_pdf_cond_mog(data, mu, Sigma, mixmat, unit_norm) +% +% Notation: Y is observation, M is mixture component, and both may be conditioned on Q. +% If Q does not exist, ignore references to Q=j below. +% Alternatively, you may ignore M if this is a conditional Gaussian. +% +% INPUTS: +% data(:,t) = t'th observation vector +% +% mu(:,k) = E[Y(t) | M(t)=k] +% or mu(:,j,k) = E[Y(t) | Q(t)=j, M(t)=k] +% +% Sigma(:,:,j,k) = Cov[Y(t) | Q(t)=j, M(t)=k] +% or there are various faster, special cases: +% Sigma() - scalar, spherical covariance independent of M,Q. +% Sigma(:,:) diag or full, tied params independent of M,Q. +% Sigma(:,:,j) tied params independent of M. +% +% mixmat(k) = Pr(M(t)=k) = prior +% or mixmat(j,k) = Pr(M(t)=k | Q(t)=j) +% Not needed if M is not defined. +% +% unit_norm - optional; if 1, means data(:,i) AND mu(:,i) each have unit norm (slightly faster) +% +% OUTPUT: +% B(t) = Pr(y(t)) +% or +% B(i,t) = Pr(y(t) | Q(t)=i) +% B2(i,k,t) = Pr(y(t) | Q(t)=i, M(t)=k) +% +% If the number of mixture components differs depending on Q, just set the trailing +% entries of mixmat to 0, e.g., 2 components if Q=1, 3 components if Q=2, +% then set mixmat(1,3)=0. In this case, B2(1,3,:)=1.0. + + + + +if isvectorBNT(mu) & size(mu,2)==1 + d = length(mu); + Q = 1; M = 1; +elseif ndims(mu)==2 + [d Q] = size(mu); + M = 1; +else + [d Q M] = size(mu); +end +[d T] = size(data); + +if nargin < 4, mixmat = ones(Q,1); end +if nargin < 5, unit_norm = 0; end + +%B2 = zeros(Q,M,T); % ATB: not needed allways +%B = zeros(Q,T); + +if isscalarBNT(Sigma) + mu = reshape(mu, [d Q*M]); + if unit_norm % (p-q)'(p-q) = p'p + q'q - 2p'q = n+m -2p'q since p(:,i)'p(:,i)=1 + %avoid an expensive repmat + disp('unit norm') + %tic; D = 2 -2*(data'*mu)'; toc + D = 2 - 2*(mu'*data); + tic; D2 = sqdist(data, mu)'; toc + assert(approxeq(D,D2)) + else + D = sqdist(data, mu)'; + end + clear mu data % ATB: clear big old data + % D(qm,t) = sq dist between data(:,t) and mu(:,qm) + logB2 = -(d/2)*log(2*pi*Sigma) - (1/(2*Sigma))*D; % det(sigma*I) = sigma^d + B2 = reshape(exp(logB2), [Q M T]); + clear logB2 % ATB: clear big old data + +elseif ndims(Sigma)==2 % tied full + mu = reshape(mu, [d Q*M]); + D = sqdist(data, mu, inv(Sigma))'; + % D(qm,t) = sq dist between data(:,t) and mu(:,qm) + logB2 = -(d/2)*log(2*pi) - 0.5*logdet(Sigma) - 0.5*D; + %denom = sqrt(det(2*pi*Sigma)); + %numer = exp(-0.5 * D); + %B2 = numer/denom; + B2 = reshape(exp(logB2), [Q M T]); + +elseif ndims(Sigma)==3 % tied across M + B2 = zeros(Q,M,T); + for j=1:Q + % D(m,t) = sq dist between data(:,t) and mu(:,j,m) + if isposdef(Sigma(:,:,j)) + D = sqdist(data, permute(mu(:,j,:), [1 3 2]), inv(Sigma(:,:,j)))'; + logB2 = -(d/2)*log(2*pi) - 0.5*logdet(Sigma(:,:,j)) - 0.5*D; + B2(j,:,:) = exp(logB2); + else + error(sprintf('mixgauss_prob: Sigma(:,:,q=%d) not psd\n', j)); + end + end + +else % general case + B2 = zeros(Q,M,T); + for j=1:Q + for k=1:M + %if mixmat(j,k) > 0 + B2(j,k,:) = gaussian_prob(data, mu(:,j,k), Sigma(:,:,j,k)); + %end + end + end +end + +% B(j,t) = sum_k B2(j,k,t) * Pr(M(t)=k | Q(t)=j) + +% The repmat is actually slower than the for-loop, because it uses too much memory +% (this is true even for small T). + +%B = squeeze(sum(B2 .* repmat(mixmat, [1 1 T]), 2)); +%B = reshape(B, [Q T]); % undo effect of squeeze in case Q = 1 + +B = zeros(Q,T); +if Q < T + for q=1:Q + %B(q,:) = mixmat(q,:) * squeeze(B2(q,:,:)); % squeeze chnages order if M=1 + B(q,:) = mixmat(q,:) * permute(B2(q,:,:), [2 3 1]); % vector * matrix sums over m + end +else + for t=1:T + B(:,t) = sum(mixmat .* B2(:,:,t), 2); % sum over m + end +end +%t=toc;fprintf('%5.3f\n', t) + +%tic +%A = squeeze(sum(B2 .* repmat(mixmat, [1 1 T]), 2)); +%t=toc;fprintf('%5.3f\n', t) +%assert(approxeq(A,B)) % may be false because of round off error diff --git a/sourcecodes/bnt-master/KPMstats/mixgauss_prob_test.m b/sourcecodes/bnt-master/KPMstats/mixgauss_prob_test.m new file mode 100644 index 00000000..a93eaf97 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/mixgauss_prob_test.m @@ -0,0 +1,111 @@ +function test_eval_pdf_cond_mixgauss() + +%Q = 10; M = 100; d = 20; T = 500; +Q = 2; M = 3; d = 4; T = 5; + +mu = rand(d,Q,M); +data = randn(d,T); +%mixmat = mk_stochastic(rand(Q,M)); +mixmat = mk_stochastic(ones(Q,M)); + +% tied scalar +Sigma = 0.01; + +mu = rand(d,M,Q); +weights = mixmat'; +N = M*ones(1,Q); +tic; [B, B2, D] = parzen(data, mu, Sigma, N, weights); toc +tic; [BC, B2C, DC] = parzenC(data, mu, Sigma, N); toc +approxeq(B,BC) +B2C = reshape(B2C,[M Q T]); +approxeq(B2,B2C) +DC = reshape(DC,[M Q T]); +approxeq(D,DC) + + +return + +tic; [B, B2] = eval_pdf_cond_mixgauss(data, mu, Sigma, mixmat); toc +tic; C = eval_pdf_cond_parzen(data, mu, Sigma); toc +approxeq(B,C) + +return; + + +mu = reshape(mu, [d Q*M]); + +data = mk_unit_norm(data); +mu = mk_unit_norm(mu); +tic; D = 2 -2*(data'*mu); toc % avoid an expensive repmat +tic; D2 = sqdist(data, mu); toc +approxeq(D,D2) + + +% D(t,m) = sq dist between data(:,t) and mu(:,m) +mu = reshape(mu, [d Q*M]); +D = dist2(data', mu'); +%denom = (2*pi)^(d/2)*sqrt(abs(det(C))); +denom = (2*pi*Sigma)^(d/2); % sqrt(det(2*pi*Sigma)) +numer = exp(-0.5/Sigma * D'); +B2 = numer / denom; +B2 = reshape(B2, [Q M T]); + +tic; B = squeeze(sum(B2 .* repmat(mixmat, [1 1 T]), 2)); toc + +tic +A = zeros(Q,T); +for q=1:Q + A(q,:) = mixmat(q,:) * squeeze(B2(q,:,:)); % sum over m +end +toc +assert(approxeq(A,B)) + +tic +A = zeros(Q,T); +for t=1:T + A(:,t) = sum(mixmat .* B2(:,:,t), 2); % sum over m +end +toc +assert(approxeq(A,B)) + + + + +mu = reshape(mu, [d Q M]); +B3 = zeros(Q,M,T); +for j=1:Q + for k=1:M + B3(j,k,:) = gaussian_prob(data, mu(:,j,k), Sigma*eye(d)); + end +end +assert(approxeq(B2, B3)) + +logB4 = -(d/2)*log(2*pi*Sigma) - (1/(2*Sigma))*D; % det(sigma*I) = sigma^d +B4 = reshape(exp(logB4), [Q M T]); +assert(approxeq(B4, B3)) + + + + +% tied cov matrix + +Sigma = rand_psd(d,d); +mu = reshape(mu, [d Q*M]); +D = sqdist(data, mu, inv(Sigma))'; +denom = sqrt(det(2*pi*Sigma)); +numer = exp(-0.5 * D); +B2 = numer / denom; +B2 = reshape(B2, [Q M T]); + +mu = reshape(mu, [d Q M]); +B3 = zeros(Q,M,T); +for j=1:Q + for k=1:M + B3(j,k,:) = gaussian_prob(data, mu(:,j,k), Sigma); + end +end +assert(approxeq(B2, B3)) + +logB4 = -(d/2)*log(2*pi) - 0.5*logdet(Sigma) - 0.5*D; +B4 = reshape(exp(logB4), [Q M T]); +assert(approxeq(B4, B3)) diff --git a/sourcecodes/bnt-master/KPMstats/mixgauss_sample.m b/sourcecodes/bnt-master/KPMstats/mixgauss_sample.m new file mode 100644 index 00000000..9438d5f0 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/mixgauss_sample.m @@ -0,0 +1,22 @@ +function [data, indices] = mixgauss_sample(mu, Sigma, mixweights, Nsamples) +% mixgauss_sample Sample data from a mixture of Gaussians +% function [data, indices] = mixgauss_sample(mu, Sigma, mixweights, Nsamples) +% +% Model is P(X) = sum_k mixweights(k) N(X; mu(:,k), Sigma(:,:,k)) or Sigma(k) for scalar +% data(:,i) is the i'th sample from P(X) +% indices(i) is the component from which sample i was drawn + +[D K] = size(mu); +data = zeros(D, Nsamples); +indices = sample_discrete(mixweights, 1, Nsamples); +for k=1:K + if ndims(Sigma) < 3 + sig = Sigma(k); + else + sig = Sigma(:,:,k); + end + ndx = find(indices==k); + if length(ndx) > 0 + data(:,ndx) = sample_gaussian(mu(:,k), sig, length(ndx))'; + end +end diff --git a/sourcecodes/bnt-master/KPMstats/mkPolyFvec.m b/sourcecodes/bnt-master/KPMstats/mkPolyFvec.m new file mode 100644 index 00000000..0ddbb259 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/mkPolyFvec.m @@ -0,0 +1,24 @@ +function p = mkPolyFvec(x) +% MKPOLYFVEC Make feature vector by constructing 2nd order polynomial from input data +% function p = mkPolyFvec(x) +% +% x(:,i) for example i +% p(:,i) = [x(1,i) x(2,i) x(3,i) x(1,i)^2 x(2,i)^2 x(3,i)^2 .. +% x(1,i)*x(2,i) x(1,i)*x(3,i) x(2,i)*x(3,i)]' +% +% Example +% x = [4 5 6]' +% p = [4 5 6 16 25 36 20 24 30]' + +fvec = x; +fvecSq = x.*x; +[D N] = size(x); +fvecCross = zeros(D*(D-1)/2, N); +i = 1; +for d=1:D + for d2=d+1:D + fvecCross(i,:) = x(d,:) .* x(d2,:); + i = i + 1; + end +end +p = [fvec; fvecSq; fvecCross]; diff --git a/sourcecodes/bnt-master/KPMstats/mk_unit_norm.m b/sourcecodes/bnt-master/KPMstats/mk_unit_norm.m new file mode 100644 index 00000000..99c150a1 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/mk_unit_norm.m @@ -0,0 +1,13 @@ +function B = mk_unit_norm(A) +% MK_UNIT_NORM Make each column be a unit norm vector +% function B = mk_unit_norm(A) +% +% We divide each column by its magnitude + + +[nrows ncols] = size(A); +s = sum(A.^2); +ndx = find(s==0); +s(ndx)=1; +B = A ./ repmat(sqrt(s), [nrows 1]); + diff --git a/sourcecodes/bnt-master/KPMstats/multinomial_prob.m b/sourcecodes/bnt-master/KPMstats/multinomial_prob.m new file mode 100644 index 00000000..460ec782 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/multinomial_prob.m @@ -0,0 +1,20 @@ +function B = eval_pdf_cond_multinomial(data, obsmat) +% EVAL_PDF_COND_MULTINOMIAL Evaluate pdf of conditional multinomial +% function B = eval_pdf_cond_multinomial(data, obsmat) +% +% Notation: Y = observation (O values), Q = conditioning variable (K values) +% +% Inputs: +% data(t) = t'th observation - must be an integer in {1,2,...,K}: cannot be 0! +% obsmat(i,o) = Pr(Y(t)=o | Q(t)=i) +% +% Output: +% B(i,t) = Pr(y(t) | Q(t)=i) + +[Q O] = size(obsmat); +T = prod(size(data)); % length(data); +B = zeros(Q,T); + +for t=1:T + B(:,t) = obsmat(:, data(t)); +end diff --git a/sourcecodes/bnt-master/KPMstats/multinomial_sample.m b/sourcecodes/bnt-master/KPMstats/multinomial_sample.m new file mode 100644 index 00000000..9c8d0cda --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/multinomial_sample.m @@ -0,0 +1,22 @@ +function Y = sample_cond_multinomial(X, M) +% SAMPLE_MULTINOMIAL Sample Y(i) ~ M(X(i), :) +% function Y = sample_multinomial(X, M) +% +% X(i) = i'th sample +% M(i,j) = P(Y=j | X=i) = noisy channel model +% +% e.g., if X is a binary image, +% Y = sample_multinomial(softeye(2, 0.9), X) +% will create a noisy version of X, where bits are flipped with probability 0.1 + +if any(X(:)==0) + error('data must only contain positive integers') +end + +Y = zeros(size(X)); +for i=min(X(:)):max(X(:)) + ndx = find(X==i); + Y(ndx) = sample_discrete(M(i,:), length(ndx), 1); +end + + diff --git a/sourcecodes/bnt-master/KPMstats/multipdf.m b/sourcecodes/bnt-master/KPMstats/multipdf.m new file mode 100644 index 00000000..592d0004 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/multipdf.m @@ -0,0 +1,45 @@ +function p = multipdf(x,theta) +%MULTIPDF Multinomial probability density function. +% p = multipdf(x,theta) returns the probabilities of +% vector x, under the multinomial distribution +% with parameter vector theta. +% +% Author: David Ross + +%-------------------------------------------------------- +% Check the arguments. +%-------------------------------------------------------- +error(nargchk(2,2,nargin)); + +% make sure theta is a vector +if ndims(theta) > 2 | all(size(theta) > 1) + error('theta must be a vector'); +end + +% make sure x is of the appropriate size +if ndims(x) > 2 | any(size(x) ~= size(theta)) + error('columns of X must have same length as theta'); +end + + +%-------------------------------------------------------- +% Main... +%-------------------------------------------------------- +p = prod(theta .^ x); +p = p .* factorial(sum(x)) ./ prod(factorial_v(x)); + + +%-------------------------------------------------------- +% Function factorial_v(x): computes the factorial function +% on each element of x +%-------------------------------------------------------- +function r = factorial_v(x) + +if size(x,2) == 1 + x = x'; +end + +r = []; +for y = x + r = [r factorial(y)]; +end \ No newline at end of file diff --git a/sourcecodes/bnt-master/KPMstats/multirnd.m b/sourcecodes/bnt-master/KPMstats/multirnd.m new file mode 100644 index 00000000..ba23f6fe --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/multirnd.m @@ -0,0 +1,48 @@ +function r = multirnd(theta,k) +%MULTIRND - Random vector from multinomial distribution. +% r = multirnd(theta,k) returns a vector randomly selected +% from the multinomial distribution with parameter vector +% theta, and count k (i.e. sum(r) = k). +% +% Note: if k is unspecified, then it is assumed k=1. +% +% Author: David Ross +% + +%-------------------------------------------------------- +% Check the arguments. +%-------------------------------------------------------- +error(nargchk(1,2,nargin)); + +% make sure theta is a vector +if ndims(theta) > 2 | all(size(theta) > 1) + error('theta must be a vector'); +end + +% if theta is a row vector, convert it to a column vector +if size(theta,1) == 1 + theta = theta'; +end + +% make sure k is a scalar? + +% if the number of samples has not been provided, set +% it to one +if nargin == 1 + k = 1; +end + + +%-------------------------------------------------------- +% Main... +%-------------------------------------------------------- +n = length(theta); +theta_cdf = cumsum(theta); + +r = zeros(n,1); +random_vals = rand(k,1); + +for j = 1:k + index = min(find(random_vals(j) <= theta_cdf)); + r(index) = r(index) + 1; +end \ No newline at end of file diff --git a/sourcecodes/bnt-master/KPMstats/normal_coef.m b/sourcecodes/bnt-master/KPMstats/normal_coef.m new file mode 100644 index 00000000..a185491e --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/normal_coef.m @@ -0,0 +1,7 @@ +function c = normal_coef (Sigma) +% NORMAL_COEF Compute the normalizing coefficient for a multivariate gaussian. +% c = normal_coef (Sigma) + +n = length(Sigma); +c = (2*pi)^(-n/2) * det(Sigma)^(-0.5); + diff --git a/sourcecodes/bnt-master/KPMstats/partial_corr_coef.m b/sourcecodes/bnt-master/KPMstats/partial_corr_coef.m new file mode 100644 index 00000000..a27d0fb9 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/partial_corr_coef.m @@ -0,0 +1,28 @@ +function [r, c] = partial_corr_coef(S, i, j, Y) +% PARTIAL_CORR_COEF Compute a partial correlation coefficient +% [r, c] = partial_corr_coef(S, i, j, Y) +% +% S is the covariance (or correlation) matrix for X, Y, Z +% where X=[i j], Y is conditioned on, and Z is marginalized out. +% Let S2 = Cov[X | Y] be the partial covariance matrix. +% Then c = S2(i,j) and r = c / sqrt( S2(i,i) * S2(j,j) ) +% + +% Example: Anderson (1984) p129 +% S = [1.0 0.8 -0.4; +% 0.8 1.0 -0.56; +% -0.4 -0.56 1.0]; +% r(1,3 | 2) = 0.0966 +% +% Example: Van de Geer (1971) p111 +%S = [1 0.453 0.322; +% 0.453 1.0 0.596; +% 0.322 0.596 1]; +% r(2,3 | 1) = 0.533 + +X = [i j]; +i2 = 1; % find_equiv_posns(i, X); +j2 = 2; % find_equiv_posns(j, X); +S2 = S(X,X) - S(X,Y)*inv(S(Y,Y))*S(Y,X); +c = S2(i2,j2); +r = c / sqrt(S2(i2,i2) * S2(j2,j2)); diff --git a/sourcecodes/bnt-master/KPMstats/parzen.m b/sourcecodes/bnt-master/KPMstats/parzen.m new file mode 100644 index 00000000..27a92783 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/parzen.m @@ -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 + + + diff --git a/sourcecodes/bnt-master/KPMstats/parzenC.c b/sourcecodes/bnt-master/KPMstats/parzenC.c new file mode 100644 index 00000000..8fb07f1d --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/parzenC.c @@ -0,0 +1,116 @@ +/* C mex version of parzen.m +[B,B2] = parzen(feat, mu, Sigma, Nproto); +*/ +#include "mex.h" +#include <stdio.h> +#include <math.h> + +#define PI 3.141592654 + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){ + int D, M, Q, T, d, m, q, t; + double *data, *mu, *SigmaPtr, *N, Sigma; + double *B, *dist, *B2, tmp; + const int* dim_mu; + double const1, const2, sum_m, sum_d, diff; + int Dt, DMq, Dm, MQt, Mq; + int dims_B2[3]; + + int ndim_mu, i, save_B2; + + data = mxGetPr(prhs[0]); + mu = mxGetPr(prhs[1]); + SigmaPtr = mxGetPr(prhs[2]); + Sigma = *SigmaPtr; + N = mxGetPr(prhs[3]); + + D = mxGetM(prhs[0]); + T = mxGetN(prhs[0]); + + ndim_mu = mxGetNumberOfDimensions(prhs[1]); + dim_mu = mxGetDimensions(prhs[1]); + D = dim_mu[0]; + M = dim_mu[1]; + /* printf("parzenC: nlhs=%d, D=%d, M=%d, T=%d\n", nlhs, D, M, T); */ + + /* If mu is mu(d,m,o,p), then [d M Q] = size(mu) in matlab sets Q=o*p, + i.e.. the size of all conditioning variabeles */ + Q = 1; + for (i = 2; i < ndim_mu; i++) { + /* printf("dim_mu[%d]=%d\n", i, dim_mu[i]); */ + Q = Q*dim_mu[i]; + } + + /* M = max(N) */ + M = -1000000; + for (i=0; i < Q; i++) { + /* printf("N[%d]=%d\n", i, (int) N[i]); */ + if (N[i] > M) { + M = (int) N[i]; + } + } + + /* printf("parzenC: nlhs=%d, D=%d, Q=%d, M=%d, T=%d\n", nlhs, D, Q, M, T); */ + + plhs[0] = mxCreateDoubleMatrix(Q,T, mxREAL); + B = mxGetPr(plhs[0]); + + if (nlhs >= 2) + save_B2 = 1; + else + save_B2 = 0; + + if (save_B2) { + /* printf("parzenC saving B2\n"); */ + /*plhs[1] = mxCreateDoubleMatrix(M*Q*T,1, mxREAL);*/ + dims_B2[0] = M; + dims_B2[1] = Q; + dims_B2[2] = T; + plhs[1] = mxCreateNumericArray(3, dims_B2, mxDOUBLE_CLASS, mxREAL); + B2 = mxGetPr(plhs[1]); + } else { + /* printf("parzenC not saving B2\n"); */ + } + /* + plhs[2] = mxCreateDoubleMatrix(M*Q*T,1, mxREAL); + dist = mxGetPr(plhs[2]); + */ + const1 = pow(2*PI*Sigma, -D/2.0); + const2 = -(1/(2*Sigma)); + + for (t=0; t < T; t++) { + /* printf("t=%d!\n",t); */ + Dt = D*t; + MQt = M*Q*t; + for (q=0; q < Q; q++) { + sum_m = 0; + DMq = D*M*q; + Mq = M*q; + + for (m=0; m < (int)N[q]; m++) { + sum_d = 0; + Dm = D*m; + for (d=0; d < D; d++) { + /* diff = data(d,t) - mu(d,m,q) */ + /*diff = data[d + D*t] - mu[d + D*m + D*M*q]; */ + diff = data[d + Dt] - mu[d + Dm + DMq]; + sum_d = sum_d + diff*diff; + } + /* dist[m,q,t] = dist[m + M*q + M*Q*t] = dist[m + Mq + MQt] = sum_d */ + tmp = const1 * exp(const2*sum_d); + sum_m = sum_m + tmp; + if (save_B2) + B2[m + Mq + MQt] = tmp; + } + + if (N[q]>0) { + B[q + Q*t] = (1.0/N[q]) * sum_m; + } else { + B[q + Q*t] = 0.0; + } + } + } +} + + + diff --git a/sourcecodes/bnt-master/KPMstats/parzenC.dll b/sourcecodes/bnt-master/KPMstats/parzenC.dll new file mode 100644 index 00000000..0700d413 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/parzenC.dll Binary files differdiff --git a/sourcecodes/bnt-master/KPMstats/parzenC_test.m b/sourcecodes/bnt-master/KPMstats/parzenC_test.m new file mode 100644 index 00000000..53a79681 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/parzenC_test.m @@ -0,0 +1,10 @@ +d = 2; M = 3; Q = 4; T = 5; Sigma = 10; +N = sample_discrete(normalize(ones(1,M)), 1, Q); +data = randn(d,T); +mu = randn(d,M,Q); + +[BM, B2M] = parzen(data, mu, Sigma, N); +[B, B2] = parzenC(data, mu, Sigma, N); + +approxeq(B,BM) +approxeq(B2,B2M) diff --git a/sourcecodes/bnt-master/KPMstats/parzen_fit_select_unif.m b/sourcecodes/bnt-master/KPMstats/parzen_fit_select_unif.m new file mode 100644 index 00000000..a9718475 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/parzen_fit_select_unif.m @@ -0,0 +1,45 @@ +function [mu, N, pick] = parzen_fit_select_unif(data, labels, max_proto, varargin) +% PARZEN_FIT_SELECT_UNIF Fit a parzen density estimator by selecting prototypes uniformly from data +% [mu, N, pick] = parzen_fit_select_unif(data, max_proto, labels, ...) +% +% We partition the data into different subsets based on the labels. +% We then choose up to max_proto columns from each subset, chosen uniformly. +% +% INPUTS +% data(:,t) +% labels(t) - should be in {1,2,..,Q} +% max_proto - max number of prototypes per partition +% +% Optional args +% partition_names{m} - for debugging +% boundary - do not choose prototypes which are within 'boundary' of the label transition +% +% OUTPUTS +% mu(:, m, q) for label q, prototype m for 1 <= m <= N(q) +% N(q) = number of prototypes for label q +% pick{q} = identity of the prototypes + +nclasses = max(labels); +[boundary, partition_names] = process_options(... + varargin, 'boundary', 0, 'partition_names', []); + +[D T] = size(data); +mu = zeros(D, 1, nclasses); % dynamically determine num prototypes (may be less than K) +mean_feat = mean(data,2); +pick = cell(1,nclasses); +for c=1:nclasses + ndx = find(labels==c); + if isempty(ndx) + %fprintf('no training images have label %d (%s)\n', c, partition_names{c}) + fprintf('no training images have label %d\n', c); + nviews = 1; + mu(:,1,c) = mean_feat; + else + foo = linspace(boundary+1, length(ndx-boundary), max_proto); + pick{c} = ndx(unique(floor(foo))); + nviews = length(pick{c}); + %fprintf('picking %d views for class %d=%s\n', nviews, c, class_names{c}); + mu(:,1:nviews,c) = data(:, pick{c}); + end + N(c) = nviews; +end diff --git a/sourcecodes/bnt-master/KPMstats/pca.m b/sourcecodes/bnt-master/KPMstats/pca.m new file mode 100644 index 00000000..4b7063d6 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/pca.m @@ -0,0 +1,42 @@ +function [PCcoeff, PCvec] = pca(data, N) +%PCA Principal Components Analysis +% +% Description +% PCCOEFF = PCA(DATA) computes the eigenvalues of the covariance +% matrix of the dataset DATA and returns them as PCCOEFF. These +% coefficients give the variance of DATA along the corresponding +% principal components. +% +% PCCOEFF = PCA(DATA, N) returns the largest N eigenvalues. +% +% [PCCOEFF, PCVEC] = PCA(DATA) returns the principal components as well +% as the coefficients. This is considerably more computationally +% demanding than just computing the eigenvalues. +% +% See also +% EIGDEC, GTMINIT, PPCA +% + +% Copyright (c) Ian T Nabney (1996-2001) + +if nargin == 1 + N = size(data, 2); +end + +if nargout == 1 + evals_only = logical(1); +else + evals_only = logical(0); +end + +if N ~= round(N) | N < 1 | N > size(data, 2) + error('Number of PCs must be integer, >0, < dim'); +end + +% Find the sorted eigenvalues of the data covariance matrix +if evals_only + PCcoeff = eigdec(cov(data), N); +else + [PCcoeff, PCvec] = eigdec(cov(data), N); +end + diff --git a/sourcecodes/bnt-master/KPMstats/rndcheck.m b/sourcecodes/bnt-master/KPMstats/rndcheck.m new file mode 100644 index 00000000..21432bb2 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/rndcheck.m @@ -0,0 +1,294 @@ +function [errorcode, rows, columns] = rndcheck(nargs,nparms,arg1,arg2,arg3,arg4,arg5) +%RNDCHECK error checks the argument list for the random number generators. + +% B.A. Jones 1-22-93 +% Copyright (c) 1993-98 by The MathWorks, Inc. +% $Revision: 1.1.1.1 $ $Date: 2005/04/26 02:29:22 $ + +sizeinfo = nargs - nparms; +errorcode = 0; + +if nparms == 3 + [r1 c1] = size(arg1); + [r2 c2] = size(arg2); + [r3 c3] = size(arg3); +end + +if nparms == 2 + [r1 c1] = size(arg1); + [r2 c2] = size(arg2); +end + +if sizeinfo == 0 + if nparms == 1 + [rows columns] = size(arg1); + end + + if nparms == 2 + scalararg1 = (prod(size(arg1)) == 1); + scalararg2 = (prod(size(arg2)) == 1); + if ~scalararg1 & ~scalararg2 + if r1 ~= r2 | c1 ~= c2 + errorcode = 1; + return; + end + end + if ~scalararg1 + [rows columns] = size(arg1); + elseif ~scalararg2 + [rows columns] = size(arg2); + else + [rows columns] = size(arg1); + end + end + + if nparms == 3 + scalararg1 = (prod(size(arg1)) == 1); + scalararg2 = (prod(size(arg2)) == 1); + scalararg3 = (prod(size(arg3)) == 1); + + if ~scalararg1 & ~scalararg2 + if r1 ~= r2 | c1 ~= c2 + errorcode = 1; + return; + end + end + + if ~scalararg1 & ~scalararg3 + if r1 ~= r3 | c1 ~= c3 + errorcode = 1; + return; + end + end + + if ~scalararg3 & ~scalararg2 + if r3 ~= r2 | c3 ~= c2 + errorcode = 1; + return; + end + end + if ~scalararg1 + [rows columns] = size(arg1); + elseif ~scalararg2 + [rows columns] = size(arg2); + else + [rows columns] = size(arg3); + end + end +end + +if sizeinfo == 1 + scalararg1 = (prod(size(arg1)) == 1); + if nparms == 1 + if prod(size(arg2)) ~= 2 + errorcode = 2; + return; + end + if ~scalararg1 & arg2 ~= size(arg1) + errorcode = 3; + return; + end + if (arg2(1) < 0 | arg2(2) < 0 | arg2(1) ~= round(arg2(1)) | arg2(2) ~= round(arg2(2))), + errorcode = 4; + return; + end + rows = arg2(1); + columns = arg2(2); + end + + if nparms == 2 + if prod(size(arg3)) ~= 2 + errorcode = 2; + return; + end + scalararg2 = (prod(size(arg2)) == 1); + if ~scalararg1 & ~scalararg2 + if r1 ~= r2 | c1 ~= c2 + errorcode = 1; + return; + end + end + if (arg3(1) < 0 | arg3(2) < 0 | arg3(1) ~= round(arg3(1)) | arg3(2) ~= round(arg3(2))), + errorcode = 4; + return; + end + if ~scalararg1 + if any(arg3 ~= size(arg1)) + errorcode = 3; + return; + end + [rows columns] = size(arg1); + elseif ~scalararg2 + if any(arg3 ~= size(arg2)) + errorcode = 3; + return; + end + [rows columns] = size(arg2); + else + rows = arg3(1); + columns = arg3(2); + end + end + + if nparms == 3 + if prod(size(arg4)) ~= 2 + errorcode = 2; + return; + end + scalararg1 = (prod(size(arg1)) == 1); + scalararg2 = (prod(size(arg2)) == 1); + scalararg3 = (prod(size(arg3)) == 1); + + if (arg4(1) < 0 | arg4(2) < 0 | arg4(1) ~= round(arg4(1)) | arg4(2) ~= round(arg4(2))), + errorcode = 4; + return; + end + + if ~scalararg1 & ~scalararg2 + if r1 ~= r2 | c1 ~= c2 + errorcode = 1; + return; + end + end + + if ~scalararg1 & ~scalararg3 + if r1 ~= r3 | c1 ~= c3 + errorcode = 1; + return; + end + end + + if ~scalararg3 & ~scalararg2 + if r3 ~= r2 | c3 ~= c2 + errorcode = 1; + return; + end + end + if ~scalararg1 + if any(arg4 ~= size(arg1)) + errorcode = 3; + return; + end + [rows columns] = size(arg1); + elseif ~scalararg2 + if any(arg4 ~= size(arg2)) + errorcode = 3; + return; + end + [rows columns] = size(arg2); + elseif ~scalararg3 + if any(arg4 ~= size(arg3)) + errorcode = 3; + return; + end + [rows columns] = size(arg3); + else + rows = arg4(1); + columns = arg4(2); + end + end +end + +if sizeinfo == 2 + if nparms == 1 + scalararg1 = (prod(size(arg1)) == 1); + if ~scalararg1 + [rows columns] = size(arg1); + if rows ~= arg2 | columns ~= arg3 + errorcode = 3; + return; + end + end + if (arg2 < 0 | arg3 < 0 | arg2 ~= round(arg2) | arg3 ~= round(arg3)), + errorcode = 4; + return; + end + rows = arg2; + columns = arg3; + end + + if nparms == 2 + scalararg1 = (prod(size(arg1)) == 1); + scalararg2 = (prod(size(arg2)) == 1); + if ~scalararg1 & ~scalararg2 + if r1 ~= r2 | c1 ~= c2 + errorcode = 1; + return; + end + end + if ~scalararg1 + [rows columns] = size(arg1); + if rows ~= arg3 | columns ~= arg4 + errorcode = 3; + return; + end + elseif ~scalararg2 + [rows columns] = size(arg2); + if rows ~= arg3 | columns ~= arg4 + errorcode = 3; + return; + end + else + if (arg3 < 0 | arg4 < 0 | arg3 ~= round(arg3) | arg4 ~= round(arg4)), + errorcode = 4; + return; + end + rows = arg3; + columns = arg4; + end + end + + if nparms == 3 + scalararg1 = (prod(size(arg1)) == 1); + scalararg2 = (prod(size(arg2)) == 1); + scalararg3 = (prod(size(arg3)) == 1); + + if ~scalararg1 & ~scalararg2 + if r1 ~= r2 | c1 ~= c2 + errorcode = 1; + return; + end + end + + if ~scalararg1 & ~scalararg3 + if r1 ~= r3 | c1 ~= c3 + errorcode = 1; + return; + end + end + + if ~scalararg3 & ~scalararg2 + if r3 ~= r2 | c3 ~= c2 + errorcode = 1; + return; + end + end + + if ~scalararg1 + [rows columns] = size(arg1); + if rows ~= arg4 | columns ~= arg5 + errorcode = 3; + return; + end + elseif ~scalararg2 + [rows columns] = size(arg2); + if rows ~= arg4 | columns ~= arg5 + errorcode = 3; + return; + end + elseif ~scalararg3 + [rows columns] = size(arg3); + if rows ~= arg4 | columns ~= arg5 + errorcode = 3; + return; + end + else + if (arg4 < 0 | arg5 < 0 | arg4 ~= round(arg4) | arg5 ~= round(arg5)), + errorcode = 4; + return; + end + rows = arg4; + columns = arg5; + end + end +end diff --git a/sourcecodes/bnt-master/KPMstats/sample.m b/sourcecodes/bnt-master/KPMstats/sample.m new file mode 100644 index 00000000..fc252da2 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/sample.m @@ -0,0 +1,15 @@ +function x = sample(p, n) +% SAMPLE Sample from categorical distribution. +% Returns a row vector of integers, sampled according to the probability +% distribution p. +% Uses the stick-breaking algorithm. +% Much faster algorithms are also possible. + +if nargin < 2 + n = 1; +end + +cdf = cumsum(p(:)); +for i = 1:n + x(i) = sum(cdf < rand) + 1; +end diff --git a/sourcecodes/bnt-master/KPMstats/sample_discrete.m b/sourcecodes/bnt-master/KPMstats/sample_discrete.m new file mode 100644 index 00000000..d8f6b687 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/sample_discrete.m @@ -0,0 +1,40 @@ +function M = sample_discrete(prob, r, c) +% SAMPLE_DISCRETE Like the built in 'rand', except we draw from a non-uniform discrete distrib. +% M = sample_discrete(prob, r, c) +% +% Example: sample_discrete([0.8 0.2], 1, 10) generates a row vector of 10 random integers from {1,2}, +% where the prob. of being 1 is 0.8 and the prob of being 2 is 0.2. + +n = length(prob); + +if nargin == 1 + r = 1; c = 1; +elseif nargin == 2 + c == r; +end + +R = rand(r, c); +M = ones(r, c); +cumprob = cumsum(prob(:)); + +if n < r*c + for i = 1:n-1 + M = M + (R > cumprob(i)); + end +else + % loop over the smaller index - can be much faster if length(prob) >> r*c + cumprob2 = cumprob(1:end-1); + for i=1:r + for j=1:c + M(i,j) = sum(R(i,j) > cumprob2)+1; + end + end +end + + +% Slower, even though vectorized +%cumprob = reshape(cumsum([0 prob(1:end-1)]), [1 1 n]); +%M = sum(R(:,:,ones(n,1)) > cumprob(ones(r,1),ones(c,1),:), 3); + +% convert using a binning algorithm +%M=bindex(R,cumprob); diff --git a/sourcecodes/bnt-master/KPMstats/sample_gaussian.m b/sourcecodes/bnt-master/KPMstats/sample_gaussian.m new file mode 100644 index 00000000..dbba6887 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/sample_gaussian.m @@ -0,0 +1,19 @@ +function M = sample_gaussian(mu, Sigma, N) +% SAMPLE_GAUSSIAN Draw N random row vectors from a Gaussian distribution +% samples = sample_gaussian(mean, cov, N) + +if nargin==2 + N = 1; +end + +% If Y = CX, Var(Y) = C Var(X) C'. +% So if Var(X)=I, and we want Var(Y)=Sigma, we need to find C. s.t. Sigma = C C'. +% Since Sigma is psd, we have Sigma = U D U' = (U D^0.5) (D'^0.5 U'). + +mu = mu(:); +n=length(mu); +[U,D,V] = svd(Sigma); +M = randn(n,N); +M = (U*sqrt(D))*M + mu*ones(1,N); % transform each column +M = M'; + diff --git a/sourcecodes/bnt-master/KPMstats/standardize.m b/sourcecodes/bnt-master/KPMstats/standardize.m new file mode 100644 index 00000000..18c17ec3 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/standardize.m @@ -0,0 +1,18 @@ +function [S, mu, sigma2] = standardize(M, mu, sigma2) +% function S = standardize(M, mu, sigma2) +% Make each column of M be zero mean, std 1. +% Thus each row is scaled separately. +% +% If mu, sigma2 are omitted, they are computed from M + +M = double(M); +if nargin < 2 + mu = mean(M,2); + sigma2 = std(M,0,2); + sigma2 = sigma2 + eps*(sigma2==0); +end + +[nrows ncols] = size(M); +S = M - repmat(mu(:), [1 ncols]); +S = S ./ repmat(sigma2, [1 ncols]); + diff --git a/sourcecodes/bnt-master/KPMstats/student_t_logprob.m b/sourcecodes/bnt-master/KPMstats/student_t_logprob.m new file mode 100644 index 00000000..0aea837e --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/student_t_logprob.m @@ -0,0 +1,13 @@ +function L = log_student_pdf(X, mu, lambda, alpha) +% LOG_STUDENT_PDF Evaluate the log of the multivariate student-t distribution at a point +% L = log_student_pdf(X, mu, lambda, alpha) +% +% Each column of X is evaluated. +% See Bernardo and Smith p435. + +k = length(mu); +assert(size(X,1) == k); +[k N] = size(X); +logc = gammaln(0.5*(alpha+k)) - gammaln(0.5*alpha) - (k/2)*log(alpha*pi) + 0.5*log(det(lambda)); +middle = (1 + (1/alpha)*(X-mu)'*lambda*(X-mu)); % scalar version +L = logc - ((alpha+k)/2)*log(middle); diff --git a/sourcecodes/bnt-master/KPMstats/student_t_prob.m b/sourcecodes/bnt-master/KPMstats/student_t_prob.m new file mode 100644 index 00000000..8ec32a4c --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/student_t_prob.m @@ -0,0 +1,19 @@ +function p = student_t_pdf(X, mu, lambda, alpha) +% STUDENT_T_PDF Evaluate the multivariate student-t distribution at a point +% p = student_t_pdf(X, mu, lambda, alpha) +% +% Each column of X is evaluated. +% See Bernardo and Smith p435. + +k = length(mu); +assert(size(X,1) == k); +[k N] = size(X); +numer = gamma(0.5*(alpha+k)); +denom = gamma(0.5*alpha) * (alpha*pi)^(k/2); +c = (numer/denom) * det(lambda)^(0.5); +p = c*(1 + (1/alpha)*(X-mu)'*lambda*(X-mu))^(-(alpha+k)/2); % scalar version +%m = repmat(mu(:), 1, N); +%exponent = sum((X-m)'*lambda*(X-m), 2); % column vector +%p = c*(1 + (1/alpha)*exponent).^(-(alpha+k)/2); + +keyboard diff --git a/sourcecodes/bnt-master/KPMstats/test_dir.m b/sourcecodes/bnt-master/KPMstats/test_dir.m new file mode 100644 index 00000000..24849473 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/test_dir.m @@ -0,0 +1,19 @@ +% # of sample points +n_samples = 1000; + +p = ones(3,1)/3; + +% Low Entropy +alpha = 0.5*p; + +% High Entropy +%alpha = 10*p; + +% draw n_samples random points from the 3-d dirichlet(alpha), +% and plot the results +points = zeros(3,n_samples); +for i = 1:n_samples + points(:,i) = dirichletrnd(alpha); +end + +scatter3(points(1,:)', points(2,:)', points(3,:)', 'r', '.', 'filled'); \ No newline at end of file diff --git a/sourcecodes/bnt-master/KPMstats/unidrndKPM.m b/sourcecodes/bnt-master/KPMstats/unidrndKPM.m new file mode 100644 index 00000000..2c96820d --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/unidrndKPM.m @@ -0,0 +1,7 @@ +function R = unidrndKPM(min, max, nr, nc) + +if nargin < 3 + nr = 1; nc = 1; +end + +R = unidrnd(max-min+1, nr, nc) + (min-1); diff --git a/sourcecodes/bnt-master/KPMstats/unif_discrete_sample.m b/sourcecodes/bnt-master/KPMstats/unif_discrete_sample.m new file mode 100644 index 00000000..aad19e05 --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/unif_discrete_sample.m @@ -0,0 +1,6 @@ +function r = unif_discrete_sample(n, nrows, ncols) +% UNIF_DISCRETE_SAMPLE Generate random numbers uniformly from {1,2,..,n} +% function r = unif_discrete_sample(n, nrows, ncols) +% Same as unidrnd in the stats toolbox. + +r = ceil(n .* rand(nrows,ncols)); diff --git a/sourcecodes/bnt-master/KPMstats/weightedRegression.m b/sourcecodes/bnt-master/KPMstats/weightedRegression.m new file mode 100644 index 00000000..2d2438ee --- /dev/null +++ b/sourcecodes/bnt-master/KPMstats/weightedRegression.m @@ -0,0 +1,58 @@ +function [a, b, error] = weightedRegression(x, z, w) +% [a , b, error] = fitRegression(x, z, w); +% % Weighted scalar linear regression +% +% Find a,b to minimize +% error = sum(w * |z - (a*x + b)|^2) +% and x(i) is a scalar + +if nargin < 3, w = ones(1,length(x)); end + +w = w(:)'; +x = x(:)'; +z = z(:)'; + +W = sum(w); +Y = sum(w .* z); +YY = sum(w .* z .* z); +YTY = sum(w .* z .* z); +X = sum(w .* x); +XX = sum(w .* x .* x); +XY = sum(w .* x .* z); + +[b, a] = clg_Mstep_simple(W, Y, YY, YTY, X, XX, XY); +error = sum(w .* (z - (a*x + b)).^2 ); + +if 0 + % demo + seed = 1; + rand('state', seed); randn('state', seed); + x = -10:10; + N = length(x); + noise = randn(1,N); + aTrue = rand(1,1); + bTrue = rand(1,1); + z = aTrue*x + bTrue + noise; + + w = ones(1,N); + [a, b, err] = weightedRegression(x, z, w); + + b2=regress(z(:), [x(:) ones(N,1)]); + assert(approxeq(b,b2(2))) + assert(approxeq(a,b2(1))) + + % Make sure we go through x(15) perfectly + w(15) = 1000; + [aW, bW, errW] = weightedRegression(x, z, w); + + figure; + plot(x, z, 'ro') + hold on + plot(x, a*x+b, 'bx-') + plot(x, aW*x+bW, 'gs-') + title(sprintf('a=%5.2f, aHat=%5.2f, aWHat=%5.3f, b=%5.2f, bHat=%5.2f, bWHat=%5.3f, err=%5.3f, errW=%5.3f', ... + aTrue, a, aW, bTrue, b, bW, err, errW)) + legend('truth', 'ls', 'wls') + +end + |