diff --git a/matlab/bingham_cluster.m b/matlab/bingham_cluster.m new file mode 100644 index 0000000..8f077d3 --- /dev/null +++ b/matlab/bingham_cluster.m @@ -0,0 +1,27 @@ +function [B weights] = bingham_cluster(X, min_points) +% [B weights] = bingham_cluster(X, min_points) -- where each B(i) contains fields V, Z, F + +if nargin < 2 + min_points = 20; +end + +for i=1:100 %max no of clusters + fprintf('Bingham distrubution #: %d', i); + [B(i) outliers] = bingham_fit_mlesac(X); + % outliers is a row vector + weights(i) = size(X,1) - length(outliers); + + fprintf('no of points left are', weights(i)); + % weight of a particular bingham is no of points left - no of ouliers + if length(outliers) < min_points + break; + end + X_updater = zeros(length(outliers), size(X,2)); + % initialize the outlier data for next iteration to zeros + for j = 1:length(outliers) + X_updater(j,:) = X(outliers(j),:); + end + % get them outliers + X = X_updater; + % update the outliers for the next iteration +end diff --git a/matlab/deprecated/bingham_fit_mlesac.m b/matlab/bingham_fit_mlesac.m similarity index 51% rename from matlab/deprecated/bingham_fit_mlesac.m rename to matlab/bingham_fit_mlesac.m index 3e813da..05f1aa5 100644 --- a/matlab/deprecated/bingham_fit_mlesac.m +++ b/matlab/bingham_fit_mlesac.m @@ -1,71 +1,66 @@ -function [B outliers] = bingham_fit_mlesac(X) -% [B outliers] = bingham_fit_mlesac(X) -- where B = B.{V, Z, F} - -d = size(X,1); -n = size(X,2); - -iter = 100; -p0 = 1 / surface_area_hypersphere(d-1); % uniform density for outliers -logp0 = log(p0); - -pmax = -inf; - -for i=1:iter - - % pick d points at random from X - r = randperm(n); - r = r(1:d); - Xi = X(:,r); - - % fit a Bingham to the d points - [V Z F] = bingham_fit(Xi); - %[V Z F] = bingham_fit_scatter(Xi*Xi') - - % compute data log likelihood - logp = 0; - for j=1:n - p = bingham_pdf(X(:,j), V, Z, F); - if p > p0 - logp = logp + log(p); - else - logp = logp + logp0; - end - end - - if logp > pmax - pmax = logp; - B.V = V; - B.Z = Z; - B.F = F; - - %fprintf('*** found new best with log likelihood %f ***\n', logp); - - %figure(20); - %plot_bingham_3d(V,Z,F,X'); - %figure(21); - %plot_bingham_3d_projections(V,Z,F); - - end -end - -% find inliers/outliers and fit the Bingham to all the inliers - -L = zeros(1,n); -for j=1:n - p = bingham_pdf(X(:,j), B.V, B.Z, B.F); - if p > p0 - L(j) = 1; - else - L(j) = 0; - end -end - -inliers = find(L); -outliers = find(~L); - -[B.V B.Z B.F] = bingham_fit(X(:,inliers)); - - - - - +function [B outliers] = bingham_fit_mlesac(X) +% [B outliers] = bingham_fit_mlesac(X) -- where B = B.{V, Z, F} +X = X'; +d = size(X,1); +n = size(X,2); + +iter = 100; +p0 = 1 / surface_area_hypersphere(d-1); % uniform density for outliers +logp0 = log(p0); + +pmax = -inf; + +for i=1:iter + fprintf('mlesac iteration: %d', i); + + % pick d points at random from X -> put them into X_1, X_2, X_3, X_4 + r = randperm(n); + r = r(1:d); + for j=1:d + eval(['X_' num2str(j) '= X(:,r(j));']) + end + X_combined = [X_1 X_2 X_3 X_4]; + X_combined = X_combined'; + + % fit a Bingham to the d points + bing_X_combined = bingham_fit(X_combined); + + % compute data log likelihood + logp = 0; + for j=1:n + p = bingham_pdf(X(:,j)', bing_X_combined); + if p > p0 + logp = logp + log(p); + else + logp = logp + logp0; + end + end + + % update the threshold + if logp > pmax + pmax = logp; + end + +end + +% find inliers/outliers +L = zeros(1,n); +for j=1:n + p = bingham_pdf(X(:,j)', bing_X_combined); + if p > p0 + L(j) = 1; + else + L(j) = 0; + end +end + +inliers = find(L); +outliers = find(~L); +fprintf('no of outliers were %d', outliers); + +% fit a Bingham to all the inliers +bing_return = bingham_fit(X(:,inliers)'); +B = bing_return; + + + diff --git a/matlab/deprecated/bingham_cluster.m b/matlab/deprecated/bingham_cluster.m deleted file mode 100644 index b71b759..0000000 --- a/matlab/deprecated/bingham_cluster.m +++ /dev/null @@ -1,18 +0,0 @@ -function [B cnts] = bingham_cluster(X, min_points) -% [B cnts] = bingham_cluster(X, min_points) -- where each B(i) contains fields V, Z, F - -if nargin < 2 - min_points = 20; -end - -for i=1:100 - - [B(i) outliers] = bingham_fit_mlesac(X); - cnts(i) = size(X,2) - length(outliers); - - if length(outliers) < min_points - break; - end - - X = X(:, outliers); -end