diff options
| author | ziejd2 | 2017-09-28 15:04:40 -0500 |
|---|---|---|
| committer | ziejd2 | 2017-09-28 15:04:40 -0500 |
| commit | 8070dc963753142bb86c4ed698d91fd623ed28e7 (patch) | |
| tree | d0f6dd8fc46a49b819aa55c1a90faa14d8448883 /sourcecodes/bnt-master/KPMstats | |
| parent | 7cc31810d53176e805532b2789955f4eedbce6bb (diff) | |
| download | BNW-8070dc963753142bb86c4ed698d91fd623ed28e7.tar.gz | |
BNW using Octave instead of Matlab.
This version of BNW should perform the same as the original version. The only difference is that it uses Octave instead of Matlab when running BayesNet Toolbox during parameter learning. I am calling this BNW_1.02. It can be accessed at: compbio.uthsc.edu/BNW_1.02
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 + |