first public version

.gitignore

@@ -0,0 +1,9 @@
+output
+*.mexw64
+*.mexw32
+*.mexa64
+*.mexmaci64
+
+cBinghamNorm?.m
+cBinghamGradLogNorm?.m
+cBinghamGradNormDividedByNorm?.m

examples/SE2PWNExample.m

@@ -0,0 +1,78 @@
+%% Example for using the SE2 PWN Distributioon
+
+%% Parameters
+si1squared = 0.3;
+si2squared = 1;
+si3squared = 1;
+
+rho12 = 0.9;
+rho13 = 0.8;
+rho23 = 0.5;
+si12 = rho12*sqrt(si1squared*si2squared);
+si13 = rho13*sqrt(si1squared*si3squared);
+si23 = rho23*sqrt(si2squared*si3squared);
+C = [si1squared si12 si13;
+    si12 si2squared si23;
+    si13 si23 si2squared];
+mu = [1,20,10]';
+
+%% Experiment with distribution
+% Create distribution
+p = SE2PWNDistribution(mu, C);
+
+% check covariance
+covarianceAnalytical = p.covariance4D
+covarianceNumerical = p.covariance4DNumerical
+covarianceError = covarianceAnalytical-covarianceNumerical
+
+% stocahstic sampling
+s = p.sample(50000);
+p2 = SE2PWNDistribution.fromSamples(s);
+
+%% create plot of samples
+rotateFirst = true;
+
+% apply rotation to translation
+% this is a different interpretation (first rotate, then translate)
+if rotateFirst
+    for i=1:size(s,2)
+        phi = s(1,i);
+        R = [cos(phi), -sin(phi); sin(phi), cos(phi)];
+        s(2:3, i) = R*s(2:3, i);
+    end
+end
+
+%% plot samples
+figure(1)
+quiver(s(2,:),s(3,:),cos(s(1,:)),sin(s(1,:)));
+hold on
+if rotateFirst
+    phi = mu(1);
+    R = [cos(phi), -sin(phi); sin(phi), cos(phi)];
+    muNew(1) = mu(1);
+    muNew(2:3) = R*mu(2:3);
+    quiver(muNew(2),muNew(3),cos(muNew(1)),sin(muNew(1)), 'color','red', 'linewidth', 2);
+else
+    quiver(mu(2),mu(3),cos(mu(1)),sin(mu(1)), 'color','red', 'linewidth', 2);
+end
+hold off
+axis equal
+xlabel('x')
+ylabel('y')
+
+%% plot 2D slice of density
+figure(2)
+step = 0.2;
+plotSize = 3;
+[x,y] = meshgrid([(mu(2)-plotSize):step:(mu(2)+plotSize)], [(mu(3)-plotSize):step:(mu(3)+plotSize)]);
+z = zeros(size(x));
+for i=1:size(x,1)
+    for j=1:size(x,2)
+        z(i,j) = p.pdf([mu(1), x(i,j), y(i,j)]');
+    end
+end
+surf(x,y,z)
+xlabel('x')
+ylabel('y')
+zlabel('f(\mu_1, x, y)');
+

examples/analyzeMetropolisHastings.m

@@ -0,0 +1,37 @@
+% This example analyzes the results obtain from the  Metropolis-Hastings 
+% algorithm. If the burn-in and/or skipping paramters are chosen too small, 
+% a significant correlation can be introduced in the samples.
+
+wn = WNDistribution(2,1.2);
+
+%% one sample per call
+n = 1000;
+s = zeros(1,n);
+for i=1:n
+    s(i) = wn.sampleMetropolisHastings(1);
+end
+
+% plot
+figure(1)
+scatter(s(1:end-1), s(2:end));
+setupAxisCircular('x','y');
+xlabel('sample n')
+ylabel('sample n+1')
+
+% calculate circular correlation
+twd = ToroidalWDDistribution([s(1:end-1); s(2:end)]);
+twd.circularCorrelationJammalamadaka
+
+%% all samples at once
+s2 = wn.sampleMetropolisHastings(n);
+
+% plot
+figure(2)
+scatter(s2(1:end-1), s2(2:end));
+setupAxisCircular('x','y');
+xlabel('sample n')
+ylabel('sample n+1')
+
+% calculate circular correlation
+twd2 = ToroidalWDDistribution([s2(1:end-1); s2(2:end)]);
+twd2.circularCorrelationJammalamadaka

examples/besselratioPlot.m

@@ -0,0 +1,21 @@
+% This example illustrates different approximations for the ratio of Bessel
+% functions
+
+%% Calculate values
+xBessel = 0:0.01:10;
+yBesselMatlab = arrayfun(@(a) besseli(1,a,1)/besseli(0,a,1), xBessel);
+yBesselAmos = arrayfun(@(a) besselratio(0,a), xBessel);
+yBesselApproxStienne12 = arrayfun(@(a) besselratioApprox(0,a,'stienne12'), xBessel);
+yBesselApproxStienne13 = arrayfun(@(a) besselratioApprox(0,a,'stienne13'), xBessel);
+
+%% Plot results
+plot(xBessel, yBesselAmos, 'b');
+hold on
+plot(xBessel, yBesselApproxStienne12, 'r');
+plot(xBessel, yBesselApproxStienne13, 'k');
+plot(xBessel(1:20:end), yBesselMatlab(1:20:end), 'kx');
+hold off
+xlabel('\kappa')
+ylabel('I_1(\kappa)/I_0(\kappa)');
+legend('Amos','Stienne12','Stienne13','Matlab','location','southeast')
+

examples/discreteFilter.m

@@ -0,0 +1,20 @@
+% Example for using the discrete filter
+
+% Initialize with 100 sampples
+df = DiscreteFilter(100);
+
+% Perform update
+df.updateNonlinear(@(x) exp(-(x-3)^2));
+clf
+df.wd.plot2d('b');
+hold on
+
+% Perform prediction
+df.predictNonlinear(@(x) x+0.1, WNDistribution(0,0.2))
+df.wd.plot2d('gx');
+hold off
+
+% Some plot settigns
+setupAxisCircular('x');
+legend('updated density', 'predicted density');
+xlabel('x'); ylabel('f(x)');

examples/mvnpdfbench.m

@@ -0,0 +1,104 @@
+% This function performs a benchmark of the mvnpdf and the mvnpdffast
+% functions.
+
+function mvnpdfbench
+    figure(1)
+    benchmark();
+    figure(2)
+    benchmarkSingle();
+
+    %bench(1) ;
+    %benchSingle(1);
+end
+
+function benchmark
+    dims = 1:20;
+    evals = 1000000;
+    repeats = 10;
+    t1 = zeros(1,repeats);
+    t2 = zeros(1,repeats);
+    times = zeros(2, length(dims));
+    for i = dims;
+        for j=1:repeats
+            [t1(j), t2(j)]= bench(i, evals);
+        end
+        times (1,i) = median(t1);
+        times (2,i) = median(t2);
+    end
+    benchPlot(dims, times, sprintf('time for %i evaluations with one call', evals))
+end
+
+function benchmarkSingle
+    dims = 1:20;
+    evals = 1000;
+    repeats = 10;
+    t1 = zeros(1,repeats);
+    t2 = zeros(1,repeats);    
+    times = zeros(2, length(dims));
+    for i = dims;
+        for j=1:repeats
+            [t1(j), t2(j)]= benchSingle(i, evals);
+        end
+        times (1,i) = median(t1);
+        times (2,i) = median(t2);
+    end
+    benchPlot(dims, times, sprintf('time for %i evaluations with individual calls', evals))
+end
+
+function benchPlot(dims, times, titletext)
+    subplot(2,1,1);
+    plot(dims, times(1,:));
+    hold on
+    plot(dims, times(2,:));
+    hold off
+    legend('mvnpdf', 'mvnpdffast', 'location', 'northwest')
+    xlabel('dimension');
+    ylabel('time (s)');
+    title(titletext);
+
+    subplot(2,1,2);
+    plot(dims, times(1,:)./times(2,:));
+    hold on
+    plot(dims, times(1,:)./times(1,:), 'b--');
+    hold off
+    xlabel('dimension');    
+    ylabel('speedup');
+end
+
+function [t1, t2]= bench(n, evals)
+    % n-D
+    rng default
+    mu = rand(1,n);
+    C = rand(n,n);
+    C=C*C';
+    x = rand(evals,n);
+
+    tic
+    mvnpdf(x, mu, C);
+    t1 = toc;
+    tic
+    mvnpdffast(x, mu, C);
+    t2 = toc;
+    fprintf('%fs\n%fs\n', t1, t2);
+end
+
+function [t1, t2]= benchSingle(n,evals)
+    % n-D
+    rng default
+    mu = rand(1,n);
+    C = rand(n,n);
+    C=C*C';
+    x = rand(evals,n);
+
+    tic
+    for i=1:length(x)
+        mvnpdf(x, mu, C);
+    end
+    t1 = toc;
+    tic
+    for i=1:length(x)
+        mvnpdffast(x, mu, C);
+    end
+    t2 = toc;
+    fprintf('%fs\n%fs\n', t1, t2);
+end

examples/stochasticSampling.m

@@ -0,0 +1,60 @@
+% This example evaluates stochastics sampling of circular densities.
+%
+% Three different sampling methods are compared, a native sample for the
+% respective density, a generic Metropolis-Hastings sampler, and a sampler
+% based on inverting the cumulative distribution function.
+
+function stochasticSampling
+    dist = WNDistribution(3, 1.5);
+    %dist = VMDistribution(3, 1.5);
+    %dist = WCDistribution(3, 1.5);
+
+    nSamples = 1000;
+
+    samples = dist.sample(nSamples);
+    figure(1);
+    compare(samples, dist, 'Native Sampling');
+
+    figure(2);
+    samples = dist.sampleMetropolisHastings(nSamples);
+    compare(samples, dist, 'Metropolis-Hastings Sampling');
+
+    figure(3);
+    samples = dist.sampleCdf(nSamples);
+    compare(samples, dist, 'Sampling Based on Inversion of CDF');
+end
+
+function compare(samples, distribution, titlestr)
+    n = size(samples,2);
+    steps = 50;
+    s = zeros(1,steps);
+    firstMomentError = zeros(1,steps);
+    kuiper = zeros(1,steps);
+    for i=1:steps
+        s(i) = floor(i/steps*n);
+        wd = WDDistribution(samples(:,1:s(i)));
+        firstMomentError(i) = abs(wd.trigonometricMoment(1) - distribution.trigonometricMoment(1));
+        kuiper(i) = wd.kuiperTest(distribution);
+    end
+    [hAx,~,~] = plotyy(s, firstMomentError, s, kuiper);
+    xlabel('samples')
+    ylabel(hAx(1),'first moment error') % left y-axis
+    ylabel(hAx(2),'kuiper test') % right y-axis
+    title(titlestr);
+
+    wd = WDDistribution(samples);
+    firstMomentErrorTotal = abs(wd.trigonometricMoment(1) - distribution.trigonometricMoment(1));
+    kuiperTotal = wd.kuiperTest(distribution);
+
+    fprintf('[%s] first moment error: %f\n', titlestr, firstMomentErrorTotal);
+    fprintf('[%s] kuiper test: %f\n', titlestr, kuiperTotal);
+
+    %llRatioWN = calculateLikelihoodRatio(samples, distribution, distribution.toWN())
+    %llRatioVM = calculateLikelihoodRatio(samples, distribution, distribution.toVM())
+    %llRatioWC = calculateLikelihoodRatio(samples, distribution, distribution.toWC())
+    %llRatioUniform = calculateLikelihoodRatio(samples, distribution, CUDistribution)
+end
+
+function lr = calculateLikelihoodRatio(samples, distribution1, distribution2)
+    lr = exp(sum(distribution1.logLikelihood(samples)) - sum(distribution2.logLikelihood(samples)));
+end

examples/wnParameterEstimation.m

@@ -0,0 +1,24 @@
+% Example for comparing different parameter estimation methods for WN
+% distributions.
+
+% Initial WN distribution
+wn = WNDistribution(4,0.4);
+
+% Obtain stochastic samples
+s = wn.sample(100);
+
+% Estimate parameters using four different methods
+wnMleJensen = WNDistribution.mleJensen(s)
+
+wnMleNumerical = WNDistribution.mleNumerical(s)
+
+wd = WDDistribution(s);
+wnMoment = wd.toWN()
+
+wnUnwrappingEM = wd.toWNunwrappingEM()
+
+% Compare logLikelihood of the samples
+wnMleJensen.logLikelihood(s)
+wnMleNumerical.logLikelihood(s)
+wnMoment.logLikelihood(s)
+wnUnwrappingEM.logLikelihood(s)

lib/compileAll.m

@@ -0,0 +1,33 @@
+% This script compiles everything.
+
+clear mex %#ok<CLSCR> % make sure mex files are not open
+
+% switch to path of this file
+cd(fileparts(mfilename('fullpath')));
+
+% verify that we are in the correct folder
+[~,dirName]=fileparts(pwd);
+if ~isequal(dirName,'lib') 
+    error('This command has to be run within the lib directory of libDirectional');
+end
+
+% util
+cd util
+compileUtil
+cd ..
+
+% externals
+cd external
+
+
+cd Faddeeva_MATLAB
+Faddeeva_build
+cd ..
+
+cd mhg13
+mex mhg.c
+mex logmhg.c
+mex mhg.c
+mex mhgi.c
+cd ..
+cd ..