From 8070dc963753142bb86c4ed698d91fd623ed28e7 Mon Sep 17 00:00:00 2001 From: ziejd2 Date: Thu, 28 Sep 2017 15:04:40 -0500 Subject: 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 --- sourcecodes/bnt-master/netlabKPM/CVS/Entries | 21 +++ sourcecodes/bnt-master/netlabKPM/CVS/Repository | 1 + sourcecodes/bnt-master/netlabKPM/CVS/Root | 1 + sourcecodes/bnt-master/netlabKPM/README.txt | 10 ++ sourcecodes/bnt-master/netlabKPM/demgmm1_movie.m | 103 ++++++++++++ .../bnt-master/netlabKPM/evidence_weighted.m | 104 ++++++++++++ sourcecodes/bnt-master/netlabKPM/glmerr_weighted.m | 57 +++++++ .../bnt-master/netlabKPM/glmgrad_weighted.m | 39 +++++ .../bnt-master/netlabKPM/glmhess_weighted.m | 57 +++++++ .../bnt-master/netlabKPM/glmtrain_weighted.m | 141 ++++++++++++++++ sourcecodes/bnt-master/netlabKPM/gmm1.avi | Bin 0 -> 93696 bytes sourcecodes/bnt-master/netlabKPM/gmmem2.m | 186 +++++++++++++++++++++ .../bnt-master/netlabKPM/gmmem_multi_restart.m | 65 +++++++ sourcecodes/bnt-master/netlabKPM/kmeans_demo.m | 76 +++++++++ sourcecodes/bnt-master/netlabKPM/mlperr_weighted.m | 68 ++++++++ .../bnt-master/netlabKPM/mlpgrad_weighted.m | 37 ++++ .../bnt-master/netlabKPM/mlphdotv_weighted.m | 73 ++++++++ .../bnt-master/netlabKPM/mlphess_weighted.m | 51 ++++++ sourcecodes/bnt-master/netlabKPM/neterr_weighted.m | 29 ++++ .../bnt-master/netlabKPM/netgrad_weighted.m | 21 +++ .../bnt-master/netlabKPM/nethess_weighted.m | 29 ++++ sourcecodes/bnt-master/netlabKPM/netopt_weighted.m | 52 ++++++ sourcecodes/bnt-master/netlabKPM/process_options.m | 132 +++++++++++++++ 23 files changed, 1353 insertions(+) create mode 100644 sourcecodes/bnt-master/netlabKPM/CVS/Entries create mode 100644 sourcecodes/bnt-master/netlabKPM/CVS/Repository create mode 100644 sourcecodes/bnt-master/netlabKPM/CVS/Root create mode 100644 sourcecodes/bnt-master/netlabKPM/README.txt create mode 100644 sourcecodes/bnt-master/netlabKPM/demgmm1_movie.m create mode 100644 sourcecodes/bnt-master/netlabKPM/evidence_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/glmerr_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/glmgrad_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/glmhess_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/glmtrain_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/gmm1.avi create mode 100644 sourcecodes/bnt-master/netlabKPM/gmmem2.m create mode 100644 sourcecodes/bnt-master/netlabKPM/gmmem_multi_restart.m create mode 100644 sourcecodes/bnt-master/netlabKPM/kmeans_demo.m create mode 100644 sourcecodes/bnt-master/netlabKPM/mlperr_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/mlpgrad_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/mlphdotv_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/mlphess_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/neterr_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/netgrad_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/nethess_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/netopt_weighted.m create mode 100644 sourcecodes/bnt-master/netlabKPM/process_options.m (limited to 'sourcecodes/bnt-master/netlabKPM') diff --git a/sourcecodes/bnt-master/netlabKPM/CVS/Entries b/sourcecodes/bnt-master/netlabKPM/CVS/Entries new file mode 100644 index 00000000..b025ed99 --- /dev/null +++ b/sourcecodes/bnt-master/netlabKPM/CVS/Entries @@ -0,0 +1,21 @@ +/README.txt/1.1.1.1/Wed Apr 27 17:59:16 2005// +/demgmm1_movie.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/evidence_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/glmerr_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/glmgrad_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/glmhess_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/glmtrain_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/gmm1.avi/1.1.1.1/Wed Apr 27 17:59:16 2005// +/gmmem2.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/gmmem_multi_restart.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/kmeans_demo.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/mlperr_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/mlpgrad_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/mlphdotv_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/mlphess_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/neterr_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/netgrad_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/nethess_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/netopt_weighted.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +/process_options.m/1.1.1.1/Wed Apr 27 17:59:16 2005// +D diff --git a/sourcecodes/bnt-master/netlabKPM/CVS/Repository b/sourcecodes/bnt-master/netlabKPM/CVS/Repository new file mode 100644 index 00000000..b51409f4 --- /dev/null +++ b/sourcecodes/bnt-master/netlabKPM/CVS/Repository @@ -0,0 +1 @@ +FullBNT/netlabKPM diff --git a/sourcecodes/bnt-master/netlabKPM/CVS/Root b/sourcecodes/bnt-master/netlabKPM/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/netlabKPM/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/netlabKPM/README.txt b/sourcecodes/bnt-master/netlabKPM/README.txt new file mode 100644 index 00000000..6ac80017 --- /dev/null +++ b/sourcecodes/bnt-master/netlabKPM/README.txt @@ -0,0 +1,10 @@ +This directory contains some small modifications/additions to +netlab, http://www.ncrg.aston.ac.uk/netlab/. + +- gmmem2 allows a prior on the covariance, and lets you visualize +intermediate results of the learning +- gmmem_multi_restart is self-explanatory + +All the _weighted functions were written by Pierpaolo Brutti. +They are needed by BNT (softmax_CPD and mlp_CPD +maximize_params method). diff --git a/sourcecodes/bnt-master/netlabKPM/demgmm1_movie.m b/sourcecodes/bnt-master/netlabKPM/demgmm1_movie.m new file mode 100644 index 00000000..783471f6 --- /dev/null +++ b/sourcecodes/bnt-master/netlabKPM/demgmm1_movie.m @@ -0,0 +1,103 @@ +%DEMGMM1 Demonstrate EM for Gaussian mixtures. +% +% Description +% This script demonstrates the use of the EM algorithm to fit a mixture +% of Gaussians to a set of data using maximum likelihood. A colour +% coding scheme is used to illustrate the evaluation of the posterior +% probabilities in the E-step of the EM algorithm. +% +% See also +% DEMGMM2, DEMGMM3, DEMGMM4, GMM, GMMEM, GMMPOST +% + +% Copyright (c) Ian T Nabney (1996-2001) + +mov = avifile('movies/gmm1.avi','fps',1 ); + +% Generate the data +randn('state', 0); rand('state', 0); +gmix = gmm(2, 2, 'spherical'); +ndat1 = 20; ndat2 = 20; ndata = ndat1+ndat2; +gmix.centres = [0.3 0.3; 0.7 0.7]; +gmix.covars = [0.01 0.01]; +x = gmmsamp(gmix, ndata); + +h = figure; +hd = plot(x(:, 1), x(:, 2), '.g', 'markersize', 30); +hold on; axis([0 1 0 1]); axis square; set(gca, 'box', 'on'); +ht = text(0.5, 1.05, 'Data', 'horizontalalignment', 'center'); + + +% Set up mixture model +ncentres = 2; input_dim = 2; +mix = gmm(input_dim, ncentres, 'spherical'); + +% Initialise the mixture model +mix.centres = [0.2 0.8; 0.8, 0.2]; +mix.covars = [0.01 0.01]; + +% Plot the initial model +ncirc = 30; theta = linspace(0, 2*pi, ncirc); +xs = cos(theta); ys = sin(theta); +xvals = mix.centres(:, 1)*ones(1,ncirc) + sqrt(mix.covars')*xs; +yvals = mix.centres(:, 2)*ones(1,ncirc) + sqrt(mix.covars')*ys; +hc(1)=line(xvals(1,:), yvals(1,:), 'color', 'r'); +hc(2)=line(xvals(2,:), yvals(2,:), 'color', 'b'); +set(ht, 'string', 'Initial Configuration'); +figure(h); +mov = addframe(mov, getframe(gcf)); +mov = addframe(mov, getframe(gcf)); + +% Initial E-step. +set(ht, 'string', 'E-step'); +post = gmmpost(mix, x); +dcols = [post(:,1), zeros(ndata, 1), post(:,2)]; +delete(hd); +for i = 1 : ndata + hd(i) = plot(x(i, 1), x(i, 2), 'color', dcols(i,:), ... + 'marker', '.', 'markersize', 30); +end + +% M-step. +set(ht, 'string', 'M-step'); +options = foptions; +options(14) = 1; % A single iteration +options(1) = -1; % Switch off all messages, including warning +mix = gmmem(mix, x, options); +delete(hc); +xvals = mix.centres(:, 1)*ones(1,ncirc) + sqrt(mix.covars')*xs; +yvals = mix.centres(:, 2)*ones(1,ncirc) + sqrt(mix.covars')*ys; +hc(1)=line(xvals(1,:), yvals(1,:), 'color', 'r'); +hc(2)=line(xvals(2,:), yvals(2,:), 'color', 'b'); +figure(h); +mov = addframe(mov, getframe(gcf)); +mov = addframe(mov, getframe(gcf)); + +% Loop over EM iterations. +numiters = 9; +for n = 1 : numiters + + set(ht, 'string', 'E-step'); + post = gmmpost(mix, x); + dcols = [post(:,1), zeros(ndata, 1), post(:,2)]; + delete(hd); + for i = 1 : ndata + hd(i) = plot(x(i, 1), x(i, 2), 'color', dcols(i,:), ... + 'marker', '.', 'markersize', 30); + end + %pause(1) + + set(ht, 'string', 'M-step'); + [mix, options] = gmmem(mix, x, options); + fprintf(1, 'Cycle %4d Error %11.6f\n', n, options(8)); + delete(hc); + xvals = mix.centres(:, 1)*ones(1,ncirc) + sqrt(mix.covars')*xs; + yvals = mix.centres(:, 2)*ones(1,ncirc) + sqrt(mix.covars')*ys; + hc(1)=line(xvals(1,:), yvals(1,:), 'color', 'r'); + hc(2)=line(xvals(2,:), yvals(2,:), 'color', 'b'); + pause(1) + + mov = addframe(mov, getframe(gcf)); +end + +mov = close(mov); diff --git a/sourcecodes/bnt-master/netlabKPM/evidence_weighted.m b/sourcecodes/bnt-master/netlabKPM/evidence_weighted.m new file mode 100644 index 00000000..3e829324 --- /dev/null +++ b/sourcecodes/bnt-master/netlabKPM/evidence_weighted.m @@ -0,0 +1,104 @@ +function [net, gamma, logev] = evidence_weighted(net, x, t, eso_w, num) +%EVIDENCE Re-estimate hyperparameters using evidence approximation. +% +% Description +% [NET] = EVIDENCE(NET, X, T) re-estimates the hyperparameters ALPHA +% and BETA by applying Bayesian re-estimation formulae for NUM +% iterations. The hyperparameter ALPHA can be a simple scalar +% associated with an isotropic prior on the weights, or can be a vector +% in which each component is associated with a group of weights as +% defined by the INDEX matrix in the NET data structure. These more +% complex priors can be set up for an MLP using MLPPRIOR. Initial +% values for the iterative re-estimation are taken from the network +% data structure NET passed as an input argument, while the return +% argument NET contains the re-estimated values. +% +% [NET, GAMMA, LOGEV] = EVIDENCE(NET, X, T, NUM) allows the re- +% estimation formula to be applied for NUM cycles in which the re- +% estimated values for the hyperparameters from each cycle are used to +% re-evaluate the Hessian matrix for the next cycle. The return value +% GAMMA is the number of well-determined parameters and LOGEV is the +% log of the evidence. +% +% See also +% MLPPRIOR, NETGRAD, NETHESS, DEMEV1, DEMARD +% + +% Copyright (c) Ian T Nabney (1996-9) + +errstring = consist(net, '', x, t); +if ~isempty(errstring) + error(errstring); +end + +ndata = size(x, 1); +if nargin == 4 + num = 1; +end + +if isfield(net,'beta') + beta = net.beta; +else + beta = 1; +end; + +% Extract weights from network +pakstr = [net.type, 'pak']; +w = feval(pakstr, net); + +% Evaluate data-dependent contribution to the Hessian matrix. +[h, dh] = nethess_weighted(w, net, x, t, eso_w); + +% Now set the negative eigenvalues to zero. +[evec, evl] = eig(dh); +evl = evl.*(evl > 0); +% safe_evl is used to avoid taking log of zero +safe_evl = evl + eps.*(evl <= 0); + +% Do the re-estimation. +for k = 1 : num + [e, edata, eprior] = neterr_weighted(w, net, x, t, eso_w); + h = nethess_weighted(w, net, x, t, eso_w, dh); + % Re-estimate alpha. + if size(net.alpha) == [1 1] + % Evaluate number of well-determined parameters. + if k == 1 + % Form vector of eigenvalues + evl = diag(evl); + safe_evl = diag(safe_evl); + end + B = beta*evl; + gamma = sum(B./(B + net.alpha)); + net.alpha = 0.5*gamma/eprior; + + % Partially evaluate log evidence + logev = e - 0.5*sum(log(safe_evl)) + 0.5*net.nwts*log(net.alpha) - ... + 0.5*ndata*log(2*pi); + else + ngroups = size(net.alpha, 1); + gams = zeros(1, ngroups); + logas = zeros(1, ngroups); + traces = zeros(1, ngroups); + % Reconstruct data hessian with negative eigenvalues set to zero. + dh = evec*evl*evec'; + hinv = inv(nethess_weighted(w, net, x, t, eso_w, dh)); + for m = 1 : ngroups + group_nweights = sum(net.index(:, m)); + gams(m) = group_nweights - ... + net.alpha(m)*sum(diag(hinv).*net.index(:,m)); + net.alpha(m) = real(gams(m)/(2*eprior(m))); + % Weight alphas by number of weights in group + logas(m) = 0.5*group_nweights*log(net.alpha(m)); + % Compute sum of evalues corresponding to group + traces(m) = sum(log(safe_evl*net.index(:,m))); + end + gamma = sum(gams, 2); + logev = e - 0.5*sum(traces) + sum(logas) - 0.5*ndata*log(2*pi); + end + % Re-estimate beta. + if isfield(net, 'beta') + net.beta = 0.5*(net.nout*ndata - gamma)/edata; + end + logev = logev + 0.5*ndata*log(beta); +end + diff --git a/sourcecodes/bnt-master/netlabKPM/glmerr_weighted.m b/sourcecodes/bnt-master/netlabKPM/glmerr_weighted.m new file mode 100644 index 00000000..f0bba449 --- /dev/null +++ b/sourcecodes/bnt-master/netlabKPM/glmerr_weighted.m @@ -0,0 +1,57 @@ +function [e, edata, eprior, y, a] = glmerr_weighted(net, x, t, eso_w) +%GLMERR Evaluate error function for generalized linear model. +% +% Description +% E = GLMERR(NET, X, T) takes a generalized linear model data +% structure NET together with a matrix X of input vectors and a matrix +% T of target vectors, and evaluates the error function E. The choice +% of error function corresponds to the output unit activation function. +% Each row of X corresponds to one input vector and each row of T +% corresponds to one target vector. +% +% [E, EDATA, EPRIOR, Y, A] = GLMERR(NET, X, T) also returns the data +% and prior components of the total error. +% +% [E, EDATA, EPRIOR, Y, A] = GLMERR(NET, X) also returns a matrix Y +% giving the outputs of the models and a matrix A giving the summed +% inputs to each output unit, where each row corresponds to one +% pattern. +% +% See also +% GLM, GLMPAK, GLMUNPAK, GLMFWD, GLMGRAD, GLMTRAIN +% + +% Copyright (c) Ian T Nabney (1996-9) + +% Check arguments for consistency +errstring = consist(net, 'glm', x, t); +if ~isempty(errstring); + error(errstring); +end + +[y, a] = glmfwd(net, x); + +%switch net.actfn + switch net.outfn + + case 'softmax' % Softmax outputs + + nout = size(a,2); + % Ensure that sum(exp(a), 2) does not overflow + maxcut = log(realmax) - log(nout); + % Ensure that exp(a) > 0 + mincut = log(realmin); + a = min(a, maxcut); + a = max(a, mincut); + temp = exp(a); + y = temp./(sum(temp, 2)*ones(1,nout)); + % Ensure that log(y) is computable + y(y errold + errold = err; + w = wold; + options(8) = err; + fprintf(1, 'Error has increased: terminating\n') + return; + end + if test && n > 1 + if (max(abs(w - wold)) < options(2) && abs(err-errold) < options(3)) + options(8) = err; + return; + else + errold = err; + wold = w; + end + end +end + +options(8) = err; +if (options(1) > 0) + disp('Warning: Maximum number of iterations has been exceeded'); +end diff --git a/sourcecodes/bnt-master/netlabKPM/gmm1.avi b/sourcecodes/bnt-master/netlabKPM/gmm1.avi new file mode 100644 index 00000000..a2d66de5 Binary files /dev/null and b/sourcecodes/bnt-master/netlabKPM/gmm1.avi differ diff --git a/sourcecodes/bnt-master/netlabKPM/gmmem2.m b/sourcecodes/bnt-master/netlabKPM/gmmem2.m new file mode 100644 index 00000000..1e64036c --- /dev/null +++ b/sourcecodes/bnt-master/netlabKPM/gmmem2.m @@ -0,0 +1,186 @@ +function [mix, num_iter, ll] = gmmem_kpm(mix, x, varargin) +%GMMEM_KPM Like GMMEM, but with additional optional arguments +% function [mix, num_iter, ll] = gmmem_kpm(mix, x, varargin) +% +% Input: +% mix - structure created by gmminit or gmmem_multi_restart +% data - each row is an example +% +% Output: +% mix - modified structure +% num_iter - number of iterations needed to reach convergence +% ll - final log likelihood +% +% [ ... ] = gmmem_kpm(..., 'param1',val1, 'param2',val2, ...) allows you to +% specify optional parameter name/value pairs. +% Parameters are below [default value in brackets] +% +% 'max_iter' - maximum number of EM iterations [10] +% 'll_thresh' - change in log-likelihood threshold for convergence [1e-2] +% 'verbose' - 1 means display output while running [0] +% 'prior_cov' - this will be added to each estimated covariance +% to prevent singularities [1e-3*eye(d)] +% 'fn' - this function, if non-empty, will be called at every iteration +% (e.g., to display the parameters as they evolve) [ [] ] +% The fn is called as fn(mix, x, iter_num, fnargs). +% It is also called before the iteration starts as +% fn(mix, x, -1, fnargs), which can be used to initialize things. +% 'fnargs' - additional arguments to be passed to fn [ {} ] +% +% Modified by Kevin P Murphy, 29 Dec 2002 + + +% Check that inputs are consistent +errstring = consist(mix, 'gmm', x); +if ~isempty(errstring) + error(errstring); +end + +[ndata, xdim] = size(x); + +[max_iter, ll_thresh, verbose, prior_cov, fn, fnargs] = ... + process_options(varargin, ... + 'max_iter', 10, 'll_thresh', 1e-2, 'verbose', 1, ... + 'prior_cov', 1e-3*eye(xdim), 'fn', [], 'fnargs', {}); + +options = foptions; +if verbose, options(1)=1; else options(1)=-1; end +options(14) = max_iter; +options(3) = ll_thresh; + + +% Sort out the options +if (options(14)) + niters = options(14); +else + niters = 100; +end + +display = options(1); +test = 0; +if options(3) > 0.0 + test = 1; % Test log likelihood for termination +end + +check_covars = 0; +if options(5) >= 1 + if display >= 0 + disp('check_covars is on'); + end + check_covars = 1; % Ensure that covariances don't collapse + MIN_COVAR = eps; % Minimum singular value of covariance matrix + init_covars = mix.covars; +end + +mix0 = mix; % save init values for debugging + +if ~isempty(fn) + feval(fn, mix, x, -1, fnargs{:}); +end + +% Main loop of algorithm +for n = 1:niters + + % Calculate posteriors based on old parameters + [post, act] = gmmpost(mix, x); + + % Calculate error value if needed + if (display | test) + prob = act*(mix.priors)'; + % Error value is negative log likelihood of data + e = - sum(log(prob + eps)); + if display > 0 + fprintf(1, 'Cycle %4d Error %11.6f\n', n, e); + end + if test + if (n > 1 & abs(e - eold) < options(3)) + options(8) = e; + ll = -e; + num_iter = n; + return; %%%%%%%%%%%%%%%% Exit here if converged + else + eold = e; + end + end + end + + if ~isempty(fn) + feval(fn, mix, x, n, fnargs{:}); + end + + % Adjust the new estimates for the parameters + new_pr = sum(post, 1); + new_c = post' * x; + + % Now move new estimates to old parameter vectors + mix.priors = new_pr ./ ndata; + + mix.centres = new_c ./ (new_pr' * ones(1, mix.nin)); + + switch mix.covar_type + case 'spherical' + n2 = dist2(x, mix.centres); + for j = 1:mix.ncentres + v(j) = (post(:,j)'*n2(:,j)); + end + mix.covars = ((v./new_pr) + sum(diag(prior_cov)))./mix.nin; + if check_covars + % Ensure that no covariance is too small + for j = 1:mix.ncentres + if mix.covars(j) < MIN_COVAR + mix.covars(j) = init_covars(j); + end + end + end + case 'diag' + for j = 1:mix.ncentres + diffs = x - (ones(ndata, 1) * mix.centres(j,:)); + wts = (post(:,j)*ones(1, mix.nin)); + mix.covars(j,:) = sum((diffs.*diffs).*wts + prior_cov, 1)./new_pr(j); + end + if check_covars + % Ensure that no covariance is too small + for j = 1:mix.ncentres + if min(mix.covars(j,:)) < MIN_COVAR + mix.covars(j,:) = init_covars(j,:); + end + end + end + case 'full' + for j = 1:mix.ncentres + diffs = x - (ones(ndata, 1) * mix.centres(j,:)); + diffs = diffs.*(sqrt(post(:,j))*ones(1, mix.nin)); + mix.covars(:,:,j) = (diffs'*diffs + prior_cov)/new_pr(j); + end + if check_covars + % Ensure that no covariance is too small + for j = 1:mix.ncentres + if min(svd(mix.covars(:,:,j))) < MIN_COVAR + mix.covars(:,:,j) = init_covars(:,:,j); + end + end + end + case 'ppca' + for j = 1:mix.ncentres + diffs = x - (ones(ndata, 1) * mix.centres(j,:)); + diffs = diffs.*(sqrt(post(:,j))*ones(1, mix.nin)); + [mix.covars(j), mix.U(:,:,j), mix.lambda(j,:)] = ... + ppca((diffs'*diffs)/new_pr(j), mix.ppca_dim); + end + if check_covars + if mix.covars(j) < MIN_COVAR + mix.covars(j) = init_covars(j); + end + end + otherwise + error(['Unknown covariance type ', mix.covar_type]); + end +end + +ll = sum(log(gmmprob(mix, x))); +num_iter = n; + +%if (display >= 0) +% disp('Warning: Maximum number of iterations has been exceeded'); +%end + diff --git a/sourcecodes/bnt-master/netlabKPM/gmmem_multi_restart.m b/sourcecodes/bnt-master/netlabKPM/gmmem_multi_restart.m new file mode 100644 index 00000000..cc24b6a7 --- /dev/null +++ b/sourcecodes/bnt-master/netlabKPM/gmmem_multi_restart.m @@ -0,0 +1,65 @@ +function [means, covs, weights, ll] = gmmem_multi_restart(K, data, varargin) +% GMMEM_MULTI_RESTART Multiple restart wrapper for gmmem_kpm +% function [means, covs, weights, ll] = gmmem_multi_restart(K, data, varargin) +% +% Input: +% K = number of mixture components +% data(i,:) is the i'th example (feature vector) +% +% Output: +% The parameters for the k'th mixture component, k=1:K, are +% means(k,:), covs(:,:,k) and weights(k) +% +% [ ... ] = gmmem_multi_restart(..., 'param1',val1, 'param2',val2, ...) +% allows you to specify optional parameter name/value pairs. +% Parameters are below [default value in brackets] +% +% 'nrestarts' - number of EM restarts [2] +% 'cov_type' - 'full', 'diag' or 'spherical' ['full'] +% 'init_cov' - the initial covariance matrix [0.1*cov(data) for each k] +% 'init_means' - [] means sample from randn(); otherwise, use +% init_means(k,:,r) for the k'th comp. on the r'th restart [ [] ] +% 'restartfn' - this function, if non-empty, will be called before/after every restart +% (e.g., to display the parameters as they evolve) [ [] ] +% The fn is called as fn(mix{r}, data, restart_num, niter, outerfnargs) +% where niter is the number of iterations performed (0 initially) +% 'restartfnargs' - additional arguments to be passed to restartfn [ {} ] +% +% Optional arguments for gmmem_kpm are passed through. +% +% Written by Kevin P Murphy, 30 Dec 2002 + +[ndata nfeatures] = size(data); + +%Cinit = repmat(0.1*diag(diag(cov(data))), [1 1 K]); +Cinit = repmat(0.1*cov(data), [1 1 K]); + +[nrestarts, init_cov, init_means, cov_type, ... + restartfn, restartfnargs, unused_args] = ... + process_options(varargin, ... + 'nrestarts', 2, 'init_cov', Cinit, 'init_means', [], ... + 'cov_type', 'full', 'restartfn', [], 'restartfnargs', {}); + +mix = cell(1, nrestarts); +cost = inf*ones(1,nrestarts); + +for r=1:nrestarts + mix{r} = gmm(nfeatures, K, cov_type); % random centers + if ~isempty(init_means), mix{r}.centres = init_means(:,:,r); end + mix{r}.covars = init_cov; + if ~isempty(restartfn) + feval(restartfn, mix{r}, data, r, 0, restartfnargs{:}); + end + [mix{r}, niter, ll] = gmmem_kpm(mix{r}, data, unused_args{:}); + cost(r) = -ll; %-sum(log(gmmprob(mix{r}, data))); + if ~isempty(restartfn) + feval(restartfn, mix{r}, data, r, niter, restartfnargs{:}); + end +end + +[nll, bestr] = min(cost); +fprintf('best r = %d\n', bestr); +ll = -nll; +means = mix{bestr}.centres; +covs = mix{bestr}.covars; +weights = mix{bestr}.priors; diff --git a/sourcecodes/bnt-master/netlabKPM/kmeans_demo.m b/sourcecodes/bnt-master/netlabKPM/kmeans_demo.m new file mode 100644 index 00000000..f67d1f0a --- /dev/null +++ b/sourcecodes/bnt-master/netlabKPM/kmeans_demo.m @@ -0,0 +1,76 @@ +function kmeans_demo() + +% Generate T points from K=5 1D clusters, and try to recover the cluster +% centers using k-means. +% Requires BNT, netlab and the matlab stats toolbox v4. + +K = 5; +ndim = 1; +true_centers = 1:K; +sigma = 1e-6; +T = 100; +% data(t,:) is the t'th data point +data = zeros(T, ndim); +% ndx(t) = i means the t'th data point is sample from cluster i +%ndx = sample_discrete(normalise(ones(1,K))); +ndx = [1*ones(1,20) 2*ones(1,20) 3*ones(1,20) 4*ones(1,20) 5*ones(1,20)]; +for t=1:T + data(t) = sample_gaussian(true_centers(ndx(t)), sigma, 1); +end +plot(1:T, data, 'x') + + + +% set the centers randomly from Gauss(0) +mix = gmm(ndim, K, 'spherical'); +h = plot_centers_as_lines(mix, [], T); + +if 0 +% Place initial centers at K data points chosen at random, but add some noise +choose_ndx = randperm(T); +choose_ndx = choose_ndx(1:K); +init_centers = data(choose_ndx) + sample_gaussian(0, 0.1, K); +mix.centres = init_centers; +h = plot_centers_as_lines(mix, h, T); +end + +if 0 +% update centers using netlab k-means +options = foptions; +niter = 10; +options(14) = niter; +mix = gmminit(mix, data, options); +h = plot_centers_as_lines(mix, h, T); +end + +% use matlab stats toolbox k-means with multiple restarts +nrestarts = 5; +[idx, centers] = kmeans(data, K, 'replicates', nrestarts, ... + 'emptyAction', 'singleton', 'display', 'iter'); +mix.centres = centers; +h = plot_centers_as_lines(mix, h, T); + +% fine tune with EM; compute covariances of each cluster +options = foptions; +niter = 20; +options(1) = 1; % display cost fn at each iter +options(14) = niter; +mix = gmmem(mix, data, options); +h = plot_centers_as_lines(mix, h, T); + +%%%%%%%%% +function h = plot_centers_as_lines(mix, h, T) + +K = mix.ncentres; +hold on +if isempty(h) + for k=1:K + h(k)=line([0 T], [mix.centres(k) mix.centres(k)]); + end +else + for k=1:K + set(h(k), 'xdata', [0 T], 'ydata', [mix.centres(k) mix.centres(k)]); + end +end +hold off + diff --git a/sourcecodes/bnt-master/netlabKPM/mlperr_weighted.m b/sourcecodes/bnt-master/netlabKPM/mlperr_weighted.m new file mode 100644 index 00000000..613cc7a7 --- /dev/null +++ b/sourcecodes/bnt-master/netlabKPM/mlperr_weighted.m @@ -0,0 +1,68 @@ +function [e, edata, eprior] = mlperr_weighted(net, x, t, eso_w) +%MLPERR Evaluate error function for 2-layer network. +% +% Description +% E = MLPERR(NET, X, T) takes a network data structure NET together +% with a matrix X of input vectors and a matrix T of target vectors, +% and evaluates the error function E. The choice of error function +% corresponds to the output unit activation function. Each row of X +% corresponds to one input vector and each row of T corresponds to one +% target vector. +% +% [E, EDATA, EPRIOR] = MLPERR(NET, X, T) additionally returns the data +% and prior components of the error, assuming a zero mean Gaussian +% prior on the weights with inverse variance parameters ALPHA and BETA +% taken from the network data structure NET. +% +% See also +% MLP, MLPPAK, MLPUNPAK, MLPFWD, MLPBKP, MLPGRAD +% + +% Copyright (c) Ian T Nabney (1996-9) + +% Check arguments for consistency +errstring = consist(net, 'mlp', x, t); +if ~isempty(errstring); + error(errstring); +end +[y, z, a] = mlpfwd(net, x); + +switch net.actfn + + case 'linear' %Linear outputs + + edata = 0.5*sum(sum((y - t).^2)); + + case 'logistic' % Logistic outputs + + % Ensure that log(1-y) is computable: need exp(a) > eps + maxcut = -log(eps); + % Ensure that log(y) is computable + mincut = -log(1/realmin - 1); + a = min(a, maxcut); + a = max(a, mincut); + y = 1./(1 + exp(-a)); + edata = - sum(sum(t.*log(y) + (1 - t).*log(1 - y))); + + case 'softmax' % Softmax outputs + + nout = size(a,2); + % Ensure that sum(exp(a), 2) does not overflow + maxcut = log(realmax) - log(nout); + % Ensure that exp(a) > 0 + mincut = log(realmin); + a = min(a, maxcut); + a = max(a, mincut); + temp = exp(a); + y = temp./(sum(temp, 2)*ones(1,nout)); + % Ensure that log(y) is computable + y(y