From ba8f404e0a71ac01d5e8c1c22d6d50c2d81f83a9 Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Wed, 26 Sep 2018 16:06:38 +0200
Subject: [PATCH 01/14] refactor VBA_ExceedanceProb to include analytical
 Dirichlet

---
 utils/VBA_ExceedanceProb.m | 150 +++++++++++++++++++++++++++----------
 1 file changed, 112 insertions(+), 38 deletions(-)

diff --git a/utils/VBA_ExceedanceProb.m b/utils/VBA_ExceedanceProb.m
index f18a9082..d193713c 100644
--- a/utils/VBA_ExceedanceProb.m
+++ b/utils/VBA_ExceedanceProb.m
@@ -1,52 +1,126 @@
-function ep = VBA_ExceedanceProb(mu,Sigma,form,verbose,Nsamp)
-% calculates the exceedance probability for mutivariate Gaussian variables
-% function ep = VBA_ExceedanceProb(mu,Sigma)
+function ep = VBA_ExceedanceProb (mu, Sigma, options)
+% // VBA toolbox //////////////////////////////////////////////////////////
+%
+% ep = VBA_ExceedanceProb (mu, Sigma, options)
+%
+% Calculates the exceedance probabilities for mutivariate Gaussian or 
+% Dirichlet distributions, i.e. the probability, for each variable, to be
+% greater than all the other ones.
+%
 % IN:
 %   - mu/Sigma: sufficient statistics of the pdf
-%       -> if form='gaussian': mu=E[x] and Sigma=V[x]
-%       -> if form='dirichlet': mu=Dirichlet counts and Sigma is unused
-%   - form: 'gaussian' or 'dirichlet'
+%       -> for a Gaussian distribution, set
+%                mu = E[x]
+%                Sigma = V[x]
+%       -> for a Dirichlet distribution, set 
+%               mu = alpha, ie, the Dirichlet pseudo-counts
+%               Sigma = [], or NaN, or left undefined
+%   - options: structure with the optional fields:
+%       + verbose: display textual info {false}
+%       + method: 'sampling' or {'analytical'} derivation of the ep
+%       + nSamples: number of samples for the sampling method {1e6}
 % OUT:
-%   - ep: vector of exceedance probabilities, i.e. the probability, for
-%   each variable, to be greater than all the other ones.
+%   - ep: vector of exceedance probabilities
+%
+% Notes:
+% ~~~~~~
+%
+% The analytical solution for the Dirichlet distribution is based on 
+% numerical integration over Gamma distributions [1] derived by 
+% Joram Soch [2].
+% [1] https://arxiv.org/abs/1611.01439
+% [2] mailto:joram.soch@bccn-berlin.de
+%
+% /////////////////////////////////////////////////////////////////////////
 
-try, form; catch, form = 'gaussian'; end 
-try, verbose; catch, verbose=0; end
-try, Nsamp; catch, Nsamp=1e5; end
+% check parameters
+% =========================================================================
+
+% guess distribution
+if ~ exist ('Sigma', 'var') || isempty (Sigma) || isnan (Sigma)
+    form = 'dirichlet';
+else
+    form = 'gaussian';
+end
+
+% fill in defaults
+if ~ exist ('options', 'var')
+    options = struct;
+end
+options = VBA_check_struct (options, ...
+    'verbose', false, ...
+    'method', 'analytical', ...
+    'nSamples', 1e6 ...
+    );
+
+% initialisation
+% =========================================================================
+K = numel (mu);
+ep = ones (K, 1);
+
+% ep computation
+% =========================================================================
 
-K = size(mu,1);
-ep = ones(K,1);
-c = [1;-1];
 switch form
+    % ---------------------------------------------------------------------
     case 'gaussian'
-        r_samp = VBA_random ('Gaussian', mu, Sigma, Nsamp);
-        [~, j] = max (r_samp);
-        tmp = histc (j, 1 : length (mu));
-        ep = tmp / Nsamp;
         
-    case 'gaussian2'
-        for k=1:K
-            for l=setdiff(1:K,k)
-                ind = [k,l];
-                m = mu(ind);
-                V = Sigma(ind,ind);
-                ep(k) = ep(k)*VBA_PPM(c'*m,c'*V*c,0,0);
-            end
+        switch options.method
+            case 'sampling'
+                r_samp = VBA_random ('Gaussian', mu, Sigma, options.nSamples);
+                [~, j] = max (r_samp);
+                tmp = histc (j, 1 : length (mu));
+                ep = tmp / options.nSamples;
+        
+            case 'analytical'
+                c = [1; -1];
+                for k = 1 : K
+                    for l = setdiff (1 : K, k)
+                        ind = [k, l];
+                        m = mu(ind);
+                        V = Sigma(ind, ind);
+                        ep(k) = ep(k) * VBA_PPM (c' * m, c' * V *c, 0, 0);
+                    end
+                end
+                ep = ep ./ sum (ep);
         end
-        ep = ep./sum(ep);
         
     case 'dirichlet'
-        r_samp = VBA_random ('Dirichlet', mu, Nsamp);
-        [y, j]  = max(r_samp);
-        if any (isnan (VBA_vec (y))) % remove failed samples in limit cases
-            j(isnan (y)) = []; 
-            Nsamp = numel (j);
-            warning ('VBA_ExceedanceProb: unstable parametrization, only %d%% of samples were correctly generated.', round(100*Nsamp/numel(y)));
+    % ---------------------------------------------------------------------
+        switch options.method
+            case 'sampling'
+                r_samp = VBA_random ('Dirichlet', mu, options.nSamples);
+                [y, j]  = max(r_samp);
+                % remove failed samples in limit cases
+                if any (isnan (VBA_vec (y))) 
+                    j(isnan (y)) = []; 
+                    options.nSamples = numel (j);
+                    warning ('VBA_ExceedanceProb: unstable parametrization, only %d%% of samples were correctly generated.', round(100*Nsamp/numel(y)));
+                end
+                tmp = histc (j, 1 : length (mu));
+                ep = tmp / options.nSamples;
+                
+            case 'analytical'
+                for k = 1 : K
+                    f = @(x) integrand (x, mu(k), mu(1 : K ~= k));
+                    ep(k) = integral (f, eps, Inf);
+                end
+                ep = ep' ./ sum (ep);
         end
-        tmp = histc (j, 1 : length (mu));
-        ep = tmp / Nsamp;
-        
-    otherwise
-        error('*** VBA_ExceedanceProb: unrecognized option form = %s', form);
+     
 end
 
+end
+
+% #########################################################################
+% Integrand function for numerical integration of the Dirichlet ep
+function p = integrand (x, aj, ak)
+    p = ones(size (x));
+    % Gamma CDF
+    for k = 1 : numel (ak)
+        p = p .* gammainc (x, ak(k));                         
+    end
+    % Gamma PDF
+    p = p .* exp ((aj - 1) .* log (x) - x - gammaln (aj));         
+end               
+

From bc0252f1765190d100bb6ca9b16eca227c3f4cd5 Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Wed, 26 Sep 2018 16:16:23 +0200
Subject: [PATCH 02/14] rename call to VBA_exceedanceProb

---
 VBA_groupBMC.m | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/VBA_groupBMC.m b/VBA_groupBMC.m
index 5c7708ca..c4b08f2e 100644
--- a/VBA_groupBMC.m
+++ b/VBA_groupBMC.m
@@ -169,9 +169,9 @@
 %-- wrap up VBA output
 out = wrapUp(L,posterior,priors,F,options);
 try
-    out.ep = VBA_ExceedanceProb(posterior.a,[],'dirichlet',0);
+    out.ep = VBA_ExceedanceProb(posterior.a);
     if ~isempty(out.options.families)
-        out.families.ep = VBA_ExceedanceProb(out.families.a,[],'dirichlet',0);
+        out.families.ep = VBA_ExceedanceProb(out.families.a);
     end
 catch
     if options.verbose
@@ -218,8 +218,7 @@
 % derive first and second order moments on model frequencies:
 [out.Ef,out.Vf] = VBA_dirichlet_moments(posterior.a);
 % derive exceedance probabilities
-% out.ep = VBA_ExceedanceProb(out.Ef,out.Vf,'gaussian');
-out.ep = VBA_ExceedanceProb(posterior.a,[],'dirichlet',0);
+out.ep = VBA_ExceedanceProb(posterior.a);
 % store accuracy and entropy terms of the Free Energy
 [F,out.ELJ,out.Sqf,out.Sqm] = FE(L,posterior,priors);
 % derive Free Energy under the null:
@@ -238,8 +237,7 @@
     out.families.r = options.C'*posterior.r;
     out.families.a = options.C'*posterior.a;
     [out.families.Ef,out.families.Vf] = VBA_dirichlet_moments(out.families.a);
-%     out.families.ep = VBA_ExceedanceProb(out.families.Ef,out.families.Vf,'gaussian');
-    out.families.ep = VBA_ExceedanceProb(out.families.a,[],'dirichlet',0);
+    out.families.ep = VBA_ExceedanceProb(out.families.a);
 end
 
 

From 6041d36f2164c7c4e6957de8228668b21718cd69 Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Wed, 26 Sep 2018 16:21:05 +0200
Subject: [PATCH 03/14] normalize file name

---
 VBA_groupBMC.m                                            | 8 ++++----
 .../{VBA_ExceedanceProb.m => VBA_exceedanceProbability.m} | 6 +++---
 2 files changed, 7 insertions(+), 7 deletions(-)
 rename utils/{VBA_ExceedanceProb.m => VBA_exceedanceProbability.m} (93%)

diff --git a/VBA_groupBMC.m b/VBA_groupBMC.m
index c4b08f2e..9bbe3610 100644
--- a/VBA_groupBMC.m
+++ b/VBA_groupBMC.m
@@ -169,9 +169,9 @@
 %-- wrap up VBA output
 out = wrapUp(L,posterior,priors,F,options);
 try
-    out.ep = VBA_ExceedanceProb(posterior.a);
+    out.ep = VBA_exceedanceProbability (posterior.a);
     if ~isempty(out.options.families)
-        out.families.ep = VBA_ExceedanceProb(out.families.a);
+        out.families.ep = VBA_exceedanceProbability (out.families.a);
     end
 catch
     if options.verbose
@@ -218,7 +218,7 @@
 % derive first and second order moments on model frequencies:
 [out.Ef,out.Vf] = VBA_dirichlet_moments(posterior.a);
 % derive exceedance probabilities
-out.ep = VBA_ExceedanceProb(posterior.a);
+out.ep = VBA_exceedanceProbability (posterior.a);
 % store accuracy and entropy terms of the Free Energy
 [F,out.ELJ,out.Sqf,out.Sqm] = FE(L,posterior,priors);
 % derive Free Energy under the null:
@@ -237,7 +237,7 @@
     out.families.r = options.C'*posterior.r;
     out.families.a = options.C'*posterior.a;
     [out.families.Ef,out.families.Vf] = VBA_dirichlet_moments(out.families.a);
-    out.families.ep = VBA_ExceedanceProb(out.families.a);
+    out.families.ep = VBA_exceedanceProbability (out.families.a);
 end
 
 
diff --git a/utils/VBA_ExceedanceProb.m b/utils/VBA_exceedanceProbability.m
similarity index 93%
rename from utils/VBA_ExceedanceProb.m
rename to utils/VBA_exceedanceProbability.m
index d193713c..a5b934cb 100644
--- a/utils/VBA_ExceedanceProb.m
+++ b/utils/VBA_exceedanceProbability.m
@@ -1,7 +1,7 @@
-function ep = VBA_ExceedanceProb (mu, Sigma, options)
+function ep = VBA_exceedanceProbability (mu, Sigma, options)
 % // VBA toolbox //////////////////////////////////////////////////////////
 %
-% ep = VBA_ExceedanceProb (mu, Sigma, options)
+% ep = VBA_exceedanceProbability (mu, Sigma, options)
 %
 % Calculates the exceedance probabilities for mutivariate Gaussian or 
 % Dirichlet distributions, i.e. the probability, for each variable, to be
@@ -95,7 +95,7 @@
                 if any (isnan (VBA_vec (y))) 
                     j(isnan (y)) = []; 
                     options.nSamples = numel (j);
-                    warning ('VBA_ExceedanceProb: unstable parametrization, only %d%% of samples were correctly generated.', round(100*Nsamp/numel(y)));
+                    warning ('VBA_exceedanceProbability: unstable parametrization, only %d%% of samples were correctly generated.', round(100*Nsamp/numel(y)));
                 end
                 tmp = histc (j, 1 : length (mu));
                 ep = tmp / options.nSamples;

From 1ab1f9a6c3aef69437219c7ec281e22f5d1a02eb Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Wed, 26 Sep 2018 16:36:05 +0200
Subject: [PATCH 04/14] fallback

---
 legacy/VBA_ExceedanceProb.m | 38 +++++++++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)
 create mode 100644 legacy/VBA_ExceedanceProb.m

diff --git a/legacy/VBA_ExceedanceProb.m b/legacy/VBA_ExceedanceProb.m
new file mode 100644
index 00000000..2a5d8ad7
--- /dev/null
+++ b/legacy/VBA_ExceedanceProb.m
@@ -0,0 +1,38 @@
+function ep = VBA_ExceedanceProb (mu, Sigma, form, verbose, Nsamp)
+
+% legacy code
+s = warning ('on');
+warning ('*** The function `VBA_ExceedanceProb` is now deprecated. Please see `VBA_exceedanceProbability` for an alternative.') 
+warning (s);
+
+options = struct ();
+
+try
+    form;
+catch
+    form = 'gaussian';
+end
+
+if form(end) == '2'
+    form = form(1:end-1);
+    options.method = 'analytical';
+else
+    options.method = 'sampling';
+end
+
+switch form
+    case 'gaussian'
+        assert(numel(mu) == numel(Sigma),'Sigma must be the same size as mu for Gaussian densities');
+    case 'dirichlet'
+        assert(isempty(Sigma),'Sigma must be empty for Dirichlet densities');
+end
+
+try
+    options.verbose = verbose;
+end
+
+try
+    options.nSamples = Nsamp;
+end
+
+ep = VBA_exceedanceProbability (mu, Sigma, options);
\ No newline at end of file

From 3ea427ad7c299ebb3efdbd11f8f2692752e52df5 Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Wed, 26 Sep 2018 16:44:04 +0200
Subject: [PATCH 05/14] parameters dimension checks

---
 legacy/VBA_ExceedanceProb.m       | 2 +-
 utils/VBA_exceedanceProbability.m | 4 +++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/legacy/VBA_ExceedanceProb.m b/legacy/VBA_ExceedanceProb.m
index 2a5d8ad7..ba6d5e1d 100644
--- a/legacy/VBA_ExceedanceProb.m
+++ b/legacy/VBA_ExceedanceProb.m
@@ -22,7 +22,7 @@
 
 switch form
     case 'gaussian'
-        assert(numel(mu) == numel(Sigma),'Sigma must be the same size as mu for Gaussian densities');
+        assert(~isempty(Sigma),'Sigma cannot be empty for Gaussian densities');
     case 'dirichlet'
         assert(isempty(Sigma),'Sigma must be empty for Dirichlet densities');
 end
diff --git a/utils/VBA_exceedanceProbability.m b/utils/VBA_exceedanceProbability.m
index a5b934cb..90fc5ddb 100644
--- a/utils/VBA_exceedanceProbability.m
+++ b/utils/VBA_exceedanceProbability.m
@@ -37,10 +37,12 @@
 % =========================================================================
 
 % guess distribution
-if ~ exist ('Sigma', 'var') || isempty (Sigma) || isnan (Sigma)
+if ~ exist ('Sigma', 'var') || isempty (Sigma) || VBA_isWeird(Sigma)
     form = 'dirichlet';
 else
     form = 'gaussian';
+    assert (all (size (Sigma) == numel (mu) * [1 1]), ...
+            '*** VBA_exceedanceProbability: For Gaussian densities, Sigma must be a n x n matrix if mu has n elements');
 end
 
 % fill in defaults

From e447c0bf4eb8f88b38f74accd17b14ba3054afc3 Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Wed, 26 Sep 2018 16:52:31 +0200
Subject: [PATCH 06/14] setup test function

---
 tests/utils/test_VBA_exceedanceProbability.m | 11 +++++++++++
 1 file changed, 11 insertions(+)
 create mode 100644 tests/utils/test_VBA_exceedanceProbability.m

diff --git a/tests/utils/test_VBA_exceedanceProbability.m b/tests/utils/test_VBA_exceedanceProbability.m
new file mode 100644
index 00000000..96d3510d
--- /dev/null
+++ b/tests/utils/test_VBA_exceedanceProbability.m
@@ -0,0 +1,11 @@
+function tests = test_VBA_exceedanceProbability
+% Unit Tests for VBA_exceedanceProbability
+
+tests = functiontests (localfunctions);
+
+% -------------------------------------------------------------------------   
+function test_empty (testCase)
+    actual = VBA_exceedanceProbability ([]);
+    testCase.verifyEmpty (actual);
+    actual = VBA_exceedanceProbability ([], []);
+    testCase.verifyEmpty (actual);
\ No newline at end of file

From 09bed90ad83fd224ae16a15a9225b60db509d2e1 Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Wed, 26 Sep 2018 18:09:07 +0200
Subject: [PATCH 07/14] test for gaussian case

---
 tests/utils/test_VBA_exceedanceProbability.m | 65 +++++++++++++++++++-
 1 file changed, 62 insertions(+), 3 deletions(-)

diff --git a/tests/utils/test_VBA_exceedanceProbability.m b/tests/utils/test_VBA_exceedanceProbability.m
index 96d3510d..9e47f962 100644
--- a/tests/utils/test_VBA_exceedanceProbability.m
+++ b/tests/utils/test_VBA_exceedanceProbability.m
@@ -3,9 +3,68 @@
 
 tests = functiontests (localfunctions);
 
+% gaussian case
 % -------------------------------------------------------------------------   
-function test_empty (testCase)
+function test_gaussian_empty (testCase)
+    actual = VBA_exceedanceProbability ([], []);
+    testCase.verifyEmpty (actual);
+  
+function test_gaussian_unity (testCase)
+    % ep for univariate should be one
+    mu = rand();
+    Sigma = rand();
+    actual = VBA_exceedanceProbability (mu, Sigma);
+    testCase.verifyEqual (actual, 1);
+    
+    % ep should be the same size as mu and sum to unity
+    k = 5;
+    mu = rand(1,k);
+    Sigma = diag(randn(1,k).^2);
+    actual = VBA_exceedanceProbability (mu, Sigma);
+    testCase.verifyNumElements (actual, k);
+    testCase.verifyEqual (sum (actual), 1, 'AbsTol', eps);
+
+function test_gaussian_canonical (testCase)
+    % bivariate case, analytical solution
+    k = 2;
+    mu = randn(1, k);
+    Sigma = randn(1,k).^2;
+    actual = VBA_exceedanceProbability (mu, diag(Sigma));
+    expected = .5 * erfc (((mu(2) - mu(1)) / sqrt (2 * sum (Sigma)))); % cdf diff
+    testCase.verifyEqual (actual(1), expected, 'AbsTol', 1e-4);
+    % centered case, equal eps
+    k = 5;
+    mu = zeros(1, k);
+    Sigma = diag(randn(1,k).^2);
+    actual = VBA_exceedanceProbability (mu, Sigma);
+    expected = ones (k, 1) / k;
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-12);
+    
+ function test_gaussian_sampling (testCase)
+     options.method = 'sampling';
+     % empty
+     actual = VBA_exceedanceProbability ([], [], options);
+     testCase.verifyEmpty (actual);
+     % unity
+     actual = VBA_exceedanceProbability (1, 1, options);
+     testCase.verifyEqual (actual, 1);
+     % general case
+     k = 2;
+     mu = rand(1,k);
+     Sigma = diag(randn(1,k).^2);
+     actual = VBA_exceedanceProbability (mu, Sigma, options);
+     expected = VBA_exceedanceProbability (mu, Sigma);
+     testCase.verifyLessThan (norm (actual - expected), 1e-3);
+    
+% dirichlet case
+% -------------------------------------------------------------------------   
+function test_dirichlet_empty (testCase)
     actual = VBA_exceedanceProbability ([]);
     testCase.verifyEmpty (actual);
-    actual = VBA_exceedanceProbability ([], []);
-    testCase.verifyEmpty (actual);
\ No newline at end of file
+    
+    
+function test_dirichlet_unity (testCase)
+    mu = rand();
+    Sigma = rand();
+    actual = VBA_exceedanceProbability (mu, Sigma);
+    testCase.verifyEqual (actual, 1);
\ No newline at end of file

From e1d1268895bc2de021232d8ddce3c7a880781e62 Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Wed, 26 Sep 2018 18:09:16 +0200
Subject: [PATCH 08/14] catch trivial case

---
 utils/VBA_exceedanceProbability.m | 20 ++++++++++++++++----
 1 file changed, 16 insertions(+), 4 deletions(-)

diff --git a/utils/VBA_exceedanceProbability.m b/utils/VBA_exceedanceProbability.m
index 90fc5ddb..fe4fa327 100644
--- a/utils/VBA_exceedanceProbability.m
+++ b/utils/VBA_exceedanceProbability.m
@@ -55,6 +55,18 @@
     'nSamples', 1e6 ...
     );
 
+% catch trivial case to avoid problems with sampling
+% =========================================================================
+if isempty (mu)
+    ep = [];
+    return;
+end
+
+if isscalar (mu)
+    ep = 1;
+    return
+end
+
 % initialisation
 % =========================================================================
 K = numel (mu);
@@ -71,17 +83,17 @@
             case 'sampling'
                 r_samp = VBA_random ('Gaussian', mu, Sigma, options.nSamples);
                 [~, j] = max (r_samp);
-                tmp = histc (j, 1 : length (mu));
-                ep = tmp / options.nSamples;
+                tmp = histcounts (j, 0.5 : K + 0.5);
+                ep = tmp' / options.nSamples;
         
             case 'analytical'
                 c = [1; -1];
                 for k = 1 : K
                     for l = setdiff (1 : K, k)
                         ind = [k, l];
-                        m = mu(ind);
+                        m = VBA_vec(mu(ind));
                         V = Sigma(ind, ind);
-                        ep(k) = ep(k) * VBA_PPM (c' * m, c' * V *c, 0, 0);
+                        ep(k) = ep(k) * VBA_PPM (c' * m, c' * V *c, 0);
                     end
                 end
                 ep = ep ./ sum (ep);

From 2e012a0651b791627c50030eab2c1b3db160d013 Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Thu, 27 Sep 2018 11:33:53 +0200
Subject: [PATCH 09/14] add tests for dirichlet ep

---
 tests/utils/test_VBA_exceedanceProbability.m | 50 ++++++++++++++++++--
 utils/VBA_exceedanceProbability.m            |  8 ++--
 2 files changed, 50 insertions(+), 8 deletions(-)

diff --git a/tests/utils/test_VBA_exceedanceProbability.m b/tests/utils/test_VBA_exceedanceProbability.m
index 9e47f962..5af8d17a 100644
--- a/tests/utils/test_VBA_exceedanceProbability.m
+++ b/tests/utils/test_VBA_exceedanceProbability.m
@@ -64,7 +64,49 @@ function test_dirichlet_empty (testCase)
     
     
 function test_dirichlet_unity (testCase)
-    mu = rand();
-    Sigma = rand();
-    actual = VBA_exceedanceProbability (mu, Sigma);
-    testCase.verifyEqual (actual, 1);
\ No newline at end of file
+    alpha = rand();
+    actual = VBA_exceedanceProbability (alpha);
+    testCase.verifyEqual (actual, 1);
+    
+    % ep should be the same size as alpha and sum to unity
+    k = 5;
+    alpha = randi(5,1,k);
+    actual = VBA_exceedanceProbability (alpha);
+    testCase.verifyNumElements (actual, k);
+    testCase.verifyEqual (sum (actual), 1, 'AbsTol', eps);
+    
+function test_dirichlet_canonical (testCase)
+    % bivariate case, analytical solution
+    k = 2;
+    alpha = randi(5, 1, k);
+    actual = VBA_exceedanceProbability (alpha);
+    expected = betainc (0.5, alpha(1) , alpha(2), 'upper');
+    testCase.verifyEqual (actual(1), expected, 'AbsTol', 1e-12);
+    
+    % equifrequent case, equal eps
+    k = 5;
+    alpha = k * ones(1, k);
+    actual = VBA_exceedanceProbability (alpha);
+    expected = ones (k, 1) / k;
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-12);
+    
+    % same, with alpha < 1
+    alpha = ones(1, k)/k;
+    actual = VBA_exceedanceProbability (alpha);
+    expected = ones (k, 1) / k;
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-12);
+    
+  function test_dirichlet_sampling (testCase)
+     options.method = 'sampling';
+     % empty
+     actual = VBA_exceedanceProbability ([], NaN, options);
+     testCase.verifyEmpty (actual);
+     % unity
+     actual = VBA_exceedanceProbability (1, NaN, options);
+     testCase.verifyEqual (actual, 1);
+     % general case
+     k = 3;
+     alpha = randi(5,1,k);
+     actual = VBA_exceedanceProbability (alpha, NaN, options);
+     expected = VBA_exceedanceProbability (alpha);
+     testCase.verifyLessThan (norm (actual - expected), 1e-3);
\ No newline at end of file
diff --git a/utils/VBA_exceedanceProbability.m b/utils/VBA_exceedanceProbability.m
index fe4fa327..a7e8b3be 100644
--- a/utils/VBA_exceedanceProbability.m
+++ b/utils/VBA_exceedanceProbability.m
@@ -18,7 +18,7 @@
 %   - options: structure with the optional fields:
 %       + verbose: display textual info {false}
 %       + method: 'sampling' or {'analytical'} derivation of the ep
-%       + nSamples: number of samples for the sampling method {1e6}
+%       + nSamples: number of samples for the sampling method {5e6}
 % OUT:
 %   - ep: vector of exceedance probabilities
 %
@@ -52,7 +52,7 @@
 options = VBA_check_struct (options, ...
     'verbose', false, ...
     'method', 'analytical', ...
-    'nSamples', 1e6 ...
+    'nSamples', 5e6 ...
     );
 
 % catch trivial case to avoid problems with sampling
@@ -112,14 +112,14 @@
                     warning ('VBA_exceedanceProbability: unstable parametrization, only %d%% of samples were correctly generated.', round(100*Nsamp/numel(y)));
                 end
                 tmp = histc (j, 1 : length (mu));
-                ep = tmp / options.nSamples;
+                ep = tmp' / options.nSamples;
                 
             case 'analytical'
                 for k = 1 : K
                     f = @(x) integrand (x, mu(k), mu(1 : K ~= k));
                     ep(k) = integral (f, eps, Inf);
                 end
-                ep = ep' ./ sum (ep);
+                ep = ep ./ sum (ep);
         end
      
 end

From 3e4ad23a88666323305771906586eaf1e5d1f0da Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Thu, 27 Sep 2018 15:44:10 +0200
Subject: [PATCH 10/14] extends test + correct bug in limit case for gaussian

---
 tests/utils/test_VBA_exceedanceProbability.m | 19 ++++++++++++++++++-
 utils/VBA_exceedanceProbability.m            |  2 +-
 2 files changed, 19 insertions(+), 2 deletions(-)

diff --git a/tests/utils/test_VBA_exceedanceProbability.m b/tests/utils/test_VBA_exceedanceProbability.m
index 5af8d17a..10a5fae3 100644
--- a/tests/utils/test_VBA_exceedanceProbability.m
+++ b/tests/utils/test_VBA_exceedanceProbability.m
@@ -39,6 +39,20 @@ function test_gaussian_canonical (testCase)
     actual = VBA_exceedanceProbability (mu, Sigma);
     expected = ones (k, 1) / k;
     testCase.verifyEqual (actual, expected, 'AbsTol', 1e-12);
+    % centered case and infinite precision, equal eps
+    k = 5;
+    mu = zeros(1, k);
+    Sigma = zeros(k);
+    actual = VBA_exceedanceProbability (mu, Sigma);
+    expected = ones (k, 1) / k;
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-12);
+    % limit case, indicator ep
+    k = 5;
+    mu = [1, zeros(1, k - 1)];
+    Sigma = zeros(k);
+    actual = VBA_exceedanceProbability (mu, Sigma);
+    expected = [1 ; zeros(k - 1, 1)];
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-12);
     
  function test_gaussian_sampling (testCase)
      options.method = 'sampling';
@@ -64,6 +78,7 @@ function test_dirichlet_empty (testCase)
     
     
 function test_dirichlet_unity (testCase)
+    % ep for univariate should be one
     alpha = rand();
     actual = VBA_exceedanceProbability (alpha);
     testCase.verifyEqual (actual, 1);
@@ -96,6 +111,8 @@ function test_dirichlet_canonical (testCase)
     expected = ones (k, 1) / k;
     testCase.verifyEqual (actual, expected, 'AbsTol', 1e-12);
     
+    
+    
   function test_dirichlet_sampling (testCase)
      options.method = 'sampling';
      % empty
@@ -105,7 +122,7 @@ function test_dirichlet_sampling (testCase)
      actual = VBA_exceedanceProbability (1, NaN, options);
      testCase.verifyEqual (actual, 1);
      % general case
-     k = 3;
+     k = 5;
      alpha = randi(5,1,k);
      actual = VBA_exceedanceProbability (alpha, NaN, options);
      expected = VBA_exceedanceProbability (alpha);
diff --git a/utils/VBA_exceedanceProbability.m b/utils/VBA_exceedanceProbability.m
index a7e8b3be..b69adffe 100644
--- a/utils/VBA_exceedanceProbability.m
+++ b/utils/VBA_exceedanceProbability.m
@@ -93,7 +93,7 @@
                         ind = [k, l];
                         m = VBA_vec(mu(ind));
                         V = Sigma(ind, ind);
-                        ep(k) = ep(k) * VBA_PPM (c' * m, c' * V *c, 0);
+                        ep(k) = ep(k) * (1 - VBA_spm_Ncdf (0, c' * m, max(c' * V * c, eps)));
                     end
                 end
                 ep = ep ./ sum (ep);

From 78b8ce6b94fb038fe083a63b5374f6d3c58b2ef4 Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Thu, 27 Sep 2018 16:18:56 +0200
Subject: [PATCH 11/14] better sampling

---
 tests/utils/test_VBA_exceedanceProbability.m | 12 ++++++++---
 utils/VBA_exceedanceProbability.m            | 22 ++++++++------------
 2 files changed, 18 insertions(+), 16 deletions(-)

diff --git a/tests/utils/test_VBA_exceedanceProbability.m b/tests/utils/test_VBA_exceedanceProbability.m
index 10a5fae3..e9213d75 100644
--- a/tests/utils/test_VBA_exceedanceProbability.m
+++ b/tests/utils/test_VBA_exceedanceProbability.m
@@ -31,7 +31,7 @@ function test_gaussian_canonical (testCase)
     Sigma = randn(1,k).^2;
     actual = VBA_exceedanceProbability (mu, diag(Sigma));
     expected = .5 * erfc (((mu(2) - mu(1)) / sqrt (2 * sum (Sigma)))); % cdf diff
-    testCase.verifyEqual (actual(1), expected, 'AbsTol', 1e-4);
+    testCase.verifyEqual (actual(1), expected, 'AbsTol', 1e-12);
     % centered case, equal eps
     k = 5;
     mu = zeros(1, k);
@@ -68,7 +68,7 @@ function test_gaussian_sampling (testCase)
      Sigma = diag(randn(1,k).^2);
      actual = VBA_exceedanceProbability (mu, Sigma, options);
      expected = VBA_exceedanceProbability (mu, Sigma);
-     testCase.verifyLessThan (norm (actual - expected), 1e-3);
+     testCase.verifyEqual (actual, expected, 'AbsTol', 5e-4);
     
 % dirichlet case
 % -------------------------------------------------------------------------   
@@ -126,4 +126,10 @@ function test_dirichlet_sampling (testCase)
      alpha = randi(5,1,k);
      actual = VBA_exceedanceProbability (alpha, NaN, options);
      expected = VBA_exceedanceProbability (alpha);
-     testCase.verifyLessThan (norm (actual - expected), 1e-3);
\ No newline at end of file
+     testCase.verifyEqual (actual, expected, 'AbsTol', 5e-4);
+     % general case
+     k = 5;
+     alpha = randi(5,1,k);
+     actual = VBA_exceedanceProbability (alpha, NaN, options);
+     expected = VBA_exceedanceProbability (alpha);
+     testCase.verifyEqual (actual, expected, 'AbsTol', 5e-4);
diff --git a/utils/VBA_exceedanceProbability.m b/utils/VBA_exceedanceProbability.m
index b69adffe..0cc41dca 100644
--- a/utils/VBA_exceedanceProbability.m
+++ b/utils/VBA_exceedanceProbability.m
@@ -82,10 +82,8 @@
         switch options.method
             case 'sampling'
                 r_samp = VBA_random ('Gaussian', mu, Sigma, options.nSamples);
-                [~, j] = max (r_samp);
-                tmp = histcounts (j, 0.5 : K + 0.5);
-                ep = tmp' / options.nSamples;
-        
+                ep = sum(bsxfun (@eq, r_samp, max (r_samp)), 2);
+                ep = ep / sum (ep);
             case 'analytical'
                 c = [1; -1];
                 for k = 1 : K
@@ -104,16 +102,14 @@
         switch options.method
             case 'sampling'
                 r_samp = VBA_random ('Dirichlet', mu, options.nSamples);
-                [y, j]  = max(r_samp);
-                % remove failed samples in limit cases
-                if any (isnan (VBA_vec (y))) 
-                    j(isnan (y)) = []; 
-                    options.nSamples = numel (j);
-                    warning ('VBA_exceedanceProbability: unstable parametrization, only %d%% of samples were correctly generated.', round(100*Nsamp/numel(y)));
-                end
-                tmp = histc (j, 1 : length (mu));
-                ep = tmp' / options.nSamples;
                 
+                ep = nansum (bsxfun (@eq, r_samp, max (r_samp)), 2);
+                ep = ep / nansum (ep);
+                
+                if VBA_isWeird (r_samp) 
+                    warning ('VBA_exceedanceProbability: unstable parametrization (alpha < 1), think about trying the analytical form.');
+                end
+                 
             case 'analytical'
                 for k = 1 : K
                     f = @(x) integrand (x, mu(k), mu(1 : K ~= k));

From faf3fbf722eb0bd79133c6b85533026b63981b5d Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Fri, 28 Sep 2018 14:01:04 +0200
Subject: [PATCH 12/14] improves sampling, nalytical, and rationalise signature

---
 VBA_groupBMC.m                               |  10 +-
 legacy/VBA_ExceedanceProb.m                  |  17 +-
 tests/utils/test_VBA_exceedanceProbability.m | 214 ++++++++++++-------
 utils/VBA_exceedanceProbability.m            | 116 +++++-----
 4 files changed, 216 insertions(+), 141 deletions(-)

diff --git a/VBA_groupBMC.m b/VBA_groupBMC.m
index 9bbe3610..964c9b20 100644
--- a/VBA_groupBMC.m
+++ b/VBA_groupBMC.m
@@ -169,9 +169,9 @@
 %-- wrap up VBA output
 out = wrapUp(L,posterior,priors,F,options);
 try
-    out.ep = VBA_exceedanceProbability (posterior.a);
+    out.ep = VBA_exceedanceProbability ('Dirichlet', posterior.a);
     if ~isempty(out.options.families)
-        out.families.ep = VBA_exceedanceProbability (out.families.a);
+        out.families.ep = VBA_exceedanceProbability ('Dirichlet', out.families.a);
     end
 catch
     if options.verbose
@@ -218,7 +218,7 @@
 % derive first and second order moments on model frequencies:
 [out.Ef,out.Vf] = VBA_dirichlet_moments(posterior.a);
 % derive exceedance probabilities
-out.ep = VBA_exceedanceProbability (posterior.a);
+out.ep = VBA_exceedanceProbability ('Dirichlet', posterior.a);
 % store accuracy and entropy terms of the Free Energy
 [F,out.ELJ,out.Sqf,out.Sqm] = FE(L,posterior,priors);
 % derive Free Energy under the null:
@@ -236,8 +236,8 @@
 if ~isempty(options.families)
     out.families.r = options.C'*posterior.r;
     out.families.a = options.C'*posterior.a;
-    [out.families.Ef,out.families.Vf] = VBA_dirichlet_moments(out.families.a);
-    out.families.ep = VBA_exceedanceProbability (out.families.a);
+    [out.families.Ef,out.families.Vf] = VBA_dirichlet_moments('Dirichlet', out.families.a);
+    out.families.ep = VBA_exceedanceProbability ('Dirichlet', out.families.a);
 end
 
 
diff --git a/legacy/VBA_ExceedanceProb.m b/legacy/VBA_ExceedanceProb.m
index ba6d5e1d..5a43717c 100644
--- a/legacy/VBA_ExceedanceProb.m
+++ b/legacy/VBA_ExceedanceProb.m
@@ -20,19 +20,22 @@
     options.method = 'sampling';
 end
 
+try
+    options.verbose = verbose;
+end
+
+try
+    options.nSamples = Nsamp;
+end
+
 switch form
     case 'gaussian'
         assert(~isempty(Sigma),'Sigma cannot be empty for Gaussian densities');
+        ep = VBA_exceedanceProbability ('Gaussian', mu, Sigma, options);
     case 'dirichlet'
         assert(isempty(Sigma),'Sigma must be empty for Dirichlet densities');
+        ep = VBA_exceedanceProbability ('Dirichlet', mu, options);
 end
 
-try
-    options.verbose = verbose;
-end
 
-try
-    options.nSamples = Nsamp;
-end
 
-ep = VBA_exceedanceProbability (mu, Sigma, options);
\ No newline at end of file
diff --git a/tests/utils/test_VBA_exceedanceProbability.m b/tests/utils/test_VBA_exceedanceProbability.m
index e9213d75..21832018 100644
--- a/tests/utils/test_VBA_exceedanceProbability.m
+++ b/tests/utils/test_VBA_exceedanceProbability.m
@@ -3,133 +3,185 @@
 
 tests = functiontests (localfunctions);
 
+% general checks
+% =========================================================================
+function test_fails_on_unknown_distribution (testCase)
+    actual = @() VBA_exceedanceProbability ('ugassnai', [1 1], eye(2));
+    testCase.verifyError(actual, '')
+    
 % gaussian case
+% =========================================================================
+
+% input cheks
 % -------------------------------------------------------------------------   
-function test_gaussian_empty (testCase)
-    actual = VBA_exceedanceProbability ([], []);
-    testCase.verifyEmpty (actual);
-  
-function test_gaussian_unity (testCase)
-    % ep for univariate should be one
-    mu = rand();
-    Sigma = rand();
-    actual = VBA_exceedanceProbability (mu, Sigma);
-    testCase.verifyEqual (actual, 1);
+function test_gaussian_fails_on_empty_moment (testCase)
+    actual = @() VBA_exceedanceProbability ('Gaussian', [], []);
+    testCase.verifyError (actual, '')
     
+function test_gaussian_return_nan_for_univariate (testCase)
+    actual = VBA_exceedanceProbability ('Gaussian', 1, 1);
+    testCase.verifyEqual (actual, NaN)
+    
+% computations
+% -------------------------------------------------------------------------     
+function test_gaussian_unity_vector (testCase)
     % ep should be the same size as mu and sum to unity
     k = 5;
     mu = rand(1,k);
     Sigma = diag(randn(1,k).^2);
-    actual = VBA_exceedanceProbability (mu, Sigma);
-    testCase.verifyNumElements (actual, k);
+    actual = VBA_exceedanceProbability ('Gaussian', mu, Sigma);
+    testCase.verifySize (actual, [k, 1]);
     testCase.verifyEqual (sum (actual), 1, 'AbsTol', eps);
 
-function test_gaussian_canonical (testCase)
-    % bivariate case, analytical solution
+function test_gaussian_bivariate (testCase)
+    % bivariate case, known analöytical solution
     k = 2;
     mu = randn(1, k);
-    Sigma = randn(1,k).^2;
-    actual = VBA_exceedanceProbability (mu, diag(Sigma));
-    expected = .5 * erfc (((mu(2) - mu(1)) / sqrt (2 * sum (Sigma)))); % cdf diff
-    testCase.verifyEqual (actual(1), expected, 'AbsTol', 1e-12);
+    Sigma = diag(randn(1,k).^2);
+    % cdf of the difference of the two varaibles
+    expected = .5 * erfc (((mu(2) - mu(1)) / sqrt (2 * sum (diag( Sigma)))));
+    actual = VBA_exceedanceProbability ('Gaussian', mu, Sigma);
+    testCase.verifyEqual (actual(1), expected, 'AbsTol', 1e-11);
+    actual = VBA_exceedanceProbability ('Gaussian', mu, Sigma, struct('method','sampling'));
+    testCase.verifyEqual (actual(1), expected, 'AbsTol', 1e-3);
+    
+function test_gaussian_centered (testCase)
     % centered case, equal eps
     k = 5;
     mu = zeros(1, k);
-    Sigma = diag(randn(1,k).^2);
-    actual = VBA_exceedanceProbability (mu, Sigma);
+    Sigma = randn()^2 * eye(k);
     expected = ones (k, 1) / k;
-    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-12);
+    actual = VBA_exceedanceProbability ('Gaussian', mu, Sigma);
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-11);
+    actual = VBA_exceedanceProbability ('Gaussian', mu, Sigma, struct('method','sampling'));
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-3);
+    
+function test_gaussian_infinite_precision (testCase)
     % centered case and infinite precision, equal eps
     k = 5;
     mu = zeros(1, k);
     Sigma = zeros(k);
-    actual = VBA_exceedanceProbability (mu, Sigma);
     expected = ones (k, 1) / k;
-    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-12);
+    actual = VBA_exceedanceProbability ('Gaussian', mu, Sigma);
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-11);
+    actual = VBA_exceedanceProbability ('Gaussian', mu, Sigma, struct('method','sampling'));
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-3);
+
+function test_gaussian_limit_case (testCase)
     % limit case, indicator ep
     k = 5;
     mu = [1, zeros(1, k - 1)];
     Sigma = zeros(k);
-    actual = VBA_exceedanceProbability (mu, Sigma);
     expected = [1 ; zeros(k - 1, 1)];
-    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-12);
-    
- function test_gaussian_sampling (testCase)
-     options.method = 'sampling';
-     % empty
-     actual = VBA_exceedanceProbability ([], [], options);
-     testCase.verifyEmpty (actual);
-     % unity
-     actual = VBA_exceedanceProbability (1, 1, options);
-     testCase.verifyEqual (actual, 1);
-     % general case
-     k = 2;
-     mu = rand(1,k);
-     Sigma = diag(randn(1,k).^2);
-     actual = VBA_exceedanceProbability (mu, Sigma, options);
-     expected = VBA_exceedanceProbability (mu, Sigma);
-     testCase.verifyEqual (actual, expected, 'AbsTol', 5e-4);
+    actual = VBA_exceedanceProbability ('Gaussian', mu, Sigma);
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-11);
+    actual = VBA_exceedanceProbability ('Gaussian', mu, Sigma, struct('method','sampling'));
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-3);
+
+ function test_gaussian_sampling_approx_analytical (testCase)
+     % general case, sampling and analytical should give the same answer
+     k = 5;
+     mu = rand(1, k);
+     Sigma = diag(randn(1, k).^2);
+     actual = VBA_exceedanceProbability ('Gaussian', mu, Sigma, struct('method','sampling'));
+     expected = VBA_exceedanceProbability ('Gaussian', mu, Sigma);
+     testCase.verifyEqual (actual, expected, 'AbsTol', 1e-3);
     
 % dirichlet case
+% =========================================================================
+
+% input cheks
 % -------------------------------------------------------------------------   
-function test_dirichlet_empty (testCase)
-    actual = VBA_exceedanceProbability ([]);
-    testCase.verifyEmpty (actual);
+function test_dirichlet_fails_on_empty_moment (testCase)
+    actual = @() VBA_exceedanceProbability ('Dirichlet', []);
+    testCase.verifyError (actual, '')
     
+function test_dirichlet_return_nan_for_univariate (testCase)
+    actual = VBA_exceedanceProbability ('Dirichlet', 1);
+    testCase.verifyEqual (actual, NaN);
     
-function test_dirichlet_unity (testCase)
-    % ep for univariate should be one
-    alpha = rand();
-    actual = VBA_exceedanceProbability (alpha);
-    testCase.verifyEqual (actual, 1);
+function dirichlet_fails_on_negative (testCase)
+    actual = @() VBA_exceedanceProbability ('Dirichlet', [-1 0]);
+    testCase.verifyError (actual, '')
     
+% computations
+% -------------------------------------------------------------------------     
+function test_dirichlet_unity_vector (testCase)
     % ep should be the same size as alpha and sum to unity
     k = 5;
     alpha = randi(5,1,k);
-    actual = VBA_exceedanceProbability (alpha);
-    testCase.verifyNumElements (actual, k);
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha);
+    testCase.verifySize (actual, [k, 1]);
+    testCase.verifyEqual (sum (actual), 1, 'AbsTol', eps);
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha, struct('method','sampling'));
+    testCase.verifySize (actual, [k, 1]);
     testCase.verifyEqual (sum (actual), 1, 'AbsTol', eps);
     
-function test_dirichlet_canonical (testCase)
+function test_dirichlet_bivariate (testCase)
     % bivariate case, analytical solution
     k = 2;
     alpha = randi(5, 1, k);
-    actual = VBA_exceedanceProbability (alpha);
     expected = betainc (0.5, alpha(1) , alpha(2), 'upper');
-    testCase.verifyEqual (actual(1), expected, 'AbsTol', 1e-12);
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha);
+    testCase.verifyEqual (actual(1), expected, 'AbsTol', 1e-11);
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha, struct('method','sampling'));
+    testCase.verifyEqual (actual(1), expected, 'AbsTol', 1e-3); 
     
+function test_dirichlet_centered (testCase)
     % equifrequent case, equal eps
     k = 5;
     alpha = k * ones(1, k);
-    actual = VBA_exceedanceProbability (alpha);
     expected = ones (k, 1) / k;
-    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-12);
-    
-    % same, with alpha < 1
-    alpha = ones(1, k)/k;
-    actual = VBA_exceedanceProbability (alpha);
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha);
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-11);
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha, struct('method','sampling'));
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-3); 
+  
+function test_dirichlet_centered_small_alphas (testCase)
+    % equifrequent case, equal eps
+    k = 5;
+    alpha = ones(1, k) / k;
     expected = ones (k, 1) / k;
-    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-12);
-    
-    
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha);
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-11);
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha, struct('method','sampling'));
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-3); 
     
-  function test_dirichlet_sampling (testCase)
-     options.method = 'sampling';
-     % empty
-     actual = VBA_exceedanceProbability ([], NaN, options);
-     testCase.verifyEmpty (actual);
-     % unity
-     actual = VBA_exceedanceProbability (1, NaN, options);
-     testCase.verifyEqual (actual, 1);
-     % general case
-     k = 5;
-     alpha = randi(5,1,k);
-     actual = VBA_exceedanceProbability (alpha, NaN, options);
-     expected = VBA_exceedanceProbability (alpha);
-     testCase.verifyEqual (actual, expected, 'AbsTol', 5e-4);
-     % general case
+  function test_dirichlet_infinite_precision (testCase)
+    % centered case and infinite precision, equal eps
+    k = 5;
+    alpha = 1e4 * ones(1, k);
+    expected = ones (k, 1) / k;
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha);
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-11);
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha, struct('method','sampling'));
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-3);
+  
+  function test_dirichlet_infinite_variance (testCase)
+    % centered case and infinite precision, equal eps
+    k = 5;
+    alpha = 1e-4 * ones(1, k);
+    expected = ones (k, 1) / k;
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha);
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-11);
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha, struct('method','sampling'));
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-3); 
+
+function test_dirichlet_limit_case (testCase)
+    % limit case, indicator ep
+    k = 5;
+    alpha = [Inf, eps * ones(1, k - 1)];
+    expected = [1 ; zeros(k - 1, 1)];
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha);
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-11);
+    actual = VBA_exceedanceProbability ('Dirichlet', alpha, struct('method','sampling'));
+    testCase.verifyEqual (actual, expected, 'AbsTol', 1e-3); 
+
+ function test_dirichlet_sampling_approx_analytical (testCase)
+     % general case, sampling and analytical should give the same answer
      k = 5;
      alpha = randi(5,1,k);
-     actual = VBA_exceedanceProbability (alpha, NaN, options);
-     expected = VBA_exceedanceProbability (alpha);
-     testCase.verifyEqual (actual, expected, 'AbsTol', 5e-4);
+     actual = VBA_exceedanceProbability ('Dirichlet', alpha, struct('method','sampling'));
+     expected = VBA_exceedanceProbability ('Dirichlet', alpha);
+     testCase.verifyEqual (actual, expected, 'AbsTol', 1e-3); 
+     
\ No newline at end of file
diff --git a/utils/VBA_exceedanceProbability.m b/utils/VBA_exceedanceProbability.m
index 0cc41dca..d27ced6d 100644
--- a/utils/VBA_exceedanceProbability.m
+++ b/utils/VBA_exceedanceProbability.m
@@ -1,20 +1,25 @@
-function ep = VBA_exceedanceProbability (mu, Sigma, options)
+function ep = VBA_exceedanceProbability (distribName, varargin)
 % // VBA toolbox //////////////////////////////////////////////////////////
 %
-% ep = VBA_exceedanceProbability (mu, Sigma, options)
+% ep = VBA_exceedanceProbability (distribName, moment_1 [, moment_2] [, options])
 %
 % Calculates the exceedance probabilities for mutivariate Gaussian or 
 % Dirichlet distributions, i.e. the probability, for each variable, to be
 % greater than all the other ones.
 %
 % IN:
-%   - mu/Sigma: sufficient statistics of the pdf
+%   - distribName: type of distribution. Could be 'Gaussian' or 'Dirichlet'
+%   - moment_x: sufficient statistics of the pdf
+%
 %       -> for a Gaussian distribution, set
-%                mu = E[x]
-%                Sigma = V[x]
+%                moment_1 = E[x]
+%                moment_2 = V[x]
+%          ie, call VBA_exceedanceProbability ('Gaussian', E[x], Var[x])
+%
 %       -> for a Dirichlet distribution, set 
-%               mu = alpha, ie, the Dirichlet pseudo-counts
-%               Sigma = [], or NaN, or left undefined
+%               moment_1 = the Dirichlet pseudo-counts, alpha
+%          ie, call VBA_exceedanceProbability ('Dirichlet', alpha)
+%
 %   - options: structure with the optional fields:
 %       + verbose: display textual info {false}
 %       + method: 'sampling' or {'analytical'} derivation of the ep
@@ -36,72 +41,85 @@
 % check parameters
 % =========================================================================
 
-% guess distribution
-if ~ exist ('Sigma', 'var') || isempty (Sigma) || VBA_isWeird(Sigma)
-    form = 'dirichlet';
-else
-    form = 'gaussian';
-    assert (all (size (Sigma) == numel (mu) * [1 1]), ...
-            '*** VBA_exceedanceProbability: For Gaussian densities, Sigma must be a n x n matrix if mu has n elements');
+% check distribution and parse moments
+switch distribName
+    case 'Gaussian'
+        mu = varargin{1};
+        Sigma = varargin{2};
+        K = numel (mu);
+        assert (isvector (mu) && ismatrix (Sigma) && all (size (Sigma) == numel (mu) * [1 1]), ...
+            '*** VBA_exceedanceProbability: For Gaussian densities, mu should be a n vector and Sigma a n x n matrix.');
+    case 'Dirichlet'
+        alpha = varargin{1};
+        K = numel (alpha);
+        assert (isvector(alpha) && all (alpha > 0), ...
+            '*** VBA_exceedanceProbability: For Dirichlet densities, alpha should be a positive vector.');
+    otherwise
+        error('*** VBA_exceedanceProbability: unknown probability density type.');
 end
 
-% fill in defaults
-if ~ exist ('options', 'var')
-    options = struct;
+% check options
+parser = inputParser;
+parser.PartialMatching = false;
+
+parser.addParameter ('verbose', false, @islogical);
+parser.addParameter ('method', 'analytical', @(z) ismember(z, {'analytical','sampling'}));
+parser.addParameter ('nSamples', 1e7, @(z) z > 0 && mod(z,1) == 0);
+
+if isstruct(varargin{end})
+    parser.parse (varargin{end});
+else
+    parser.parse ();
 end
-options = VBA_check_struct (options, ...
-    'verbose', false, ...
-    'method', 'analytical', ...
-    'nSamples', 5e6 ...
-    );
+options = parser.Results ;
 
-% catch trivial case to avoid problems with sampling
+% catch trivial cases
 % =========================================================================
-if isempty (mu)
-    ep = [];
+if K < 2
+    ep = NaN;
     return;
 end
 
-if isscalar (mu)
-    ep = 1;
-    return
-end
-
 % initialisation
 % =========================================================================
-K = numel (mu);
 ep = ones (K, 1);
 
 % ep computation
 % =========================================================================
 
-switch form
+switch distribName
     % ---------------------------------------------------------------------
-    case 'gaussian'
+    case 'Gaussian'
         
         switch options.method
+            
             case 'sampling'
                 r_samp = VBA_random ('Gaussian', mu, Sigma, options.nSamples);
                 ep = sum(bsxfun (@eq, r_samp, max (r_samp)), 2);
                 ep = ep / sum (ep);
+                
             case 'analytical'
-                c = [1; -1];
+                if ~ isdiag(Sigma)
+                    error('*** VBA_exceedanceProbability: cannot compute analytical ep for non diagonal covariance. Please use the sampling method');
+                end
+                dSigma = max (diag (Sigma), 1e-6);
+                mu = VBA_vec(mu);
                 for k = 1 : K
-                    for l = setdiff (1 : K, k)
-                        ind = [k, l];
-                        m = VBA_vec(mu(ind));
-                        V = Sigma(ind, ind);
-                        ep(k) = ep(k) * (1 - VBA_spm_Ncdf (0, c' * m, max(c' * V * c, eps)));
-                    end
+                    l = setdiff (1 : K, k);  
+                    % p(x_k > x_j) = int pdf_k(t)cdf_1(t)cdf_2(t)...cdf_k-1(t)cdf_k+1(t)...cdf_k(t) dt 
+                    mycdf = @(t) sum (bsxfun(@(x,i) log(VBA_spm_Ncdf(x, mu(l(i)), dSigma(l(i)))), VBA_vec(t), 1 : (K-1)), 2)';                 
+                    ii = @(t) exp(log(VBA_spm_Npdf(t, mu(k), dSigma(k))) + mycdf(t));
+                    ep(k) = integral(ii, -Inf, Inf);
                 end
-                ep = ep ./ sum (ep);
+                ep = ep / sum(ep);
         end
-        
-    case 'dirichlet'
     % ---------------------------------------------------------------------
+    case 'Dirichlet'
+        
         switch options.method
+            
             case 'sampling'
-                r_samp = VBA_random ('Dirichlet', mu, options.nSamples);
+                r_samp = VBA_random ('Dirichlet', alpha, options.nSamples);
                 
                 ep = nansum (bsxfun (@eq, r_samp, max (r_samp)), 2);
                 ep = ep / nansum (ep);
@@ -111,11 +129,13 @@
                 end
                  
             case 'analytical'
+                alpha = min (alpha, 1e7);
+                
                 for k = 1 : K
-                    f = @(x) integrand (x, mu(k), mu(1 : K ~= k));
-                    ep(k) = integral (f, eps, Inf);
+                    f = @(x) integrand (x, alpha(k), alpha(1 : K ~= k));
+                    ep(k) = integral (f, eps, Inf, 'Waypoints', 10.^(-15:15));
                 end
-                ep = ep ./ sum (ep);
+                ep = ep / sum (ep);
         end
      
 end
@@ -131,6 +151,6 @@
         p = p .* gammainc (x, ak(k));                         
     end
     % Gamma PDF
-    p = p .* exp ((aj - 1) .* log (x) - x - gammaln (aj));         
+    p = p .* exp ((aj - 1) .* log (x) - x - gammaln (aj));    
 end               
 

From d921d791d349a1d55910a71801f76c8a0e3c889a Mon Sep 17 00:00:00 2001
From: Unknown <lionel.rigoux@sf.mpg.de>
Date: Fri, 28 Sep 2018 23:55:20 +0200
Subject: [PATCH 13/14] refine doc

---
 utils/VBA_exceedanceProbability.m | 28 +++++++++++++---------------
 1 file changed, 13 insertions(+), 15 deletions(-)

diff --git a/utils/VBA_exceedanceProbability.m b/utils/VBA_exceedanceProbability.m
index d27ced6d..f9962b47 100644
--- a/utils/VBA_exceedanceProbability.m
+++ b/utils/VBA_exceedanceProbability.m
@@ -1,7 +1,8 @@
 function ep = VBA_exceedanceProbability (distribName, varargin)
 % // VBA toolbox //////////////////////////////////////////////////////////
 %
-% ep = VBA_exceedanceProbability (distribName, moment_1 [, moment_2] [, options])
+% ep = VBA_exceedanceProbability ('Gaussian', mu, Sigma, [options])
+% ep = VBA_exceedanceProbability ('Dirichlet', alpha, [options])
 %
 % Calculates the exceedance probabilities for mutivariate Gaussian or 
 % Dirichlet distributions, i.e. the probability, for each variable, to be
@@ -9,18 +10,15 @@
 %
 % IN:
 %   - distribName: type of distribution. Could be 'Gaussian' or 'Dirichlet'
-%   - moment_x: sufficient statistics of the pdf
 %
-%       -> for a Gaussian distribution, set
-%                moment_1 = E[x]
-%                moment_2 = V[x]
-%          ie, call VBA_exceedanceProbability ('Gaussian', E[x], Var[x])
+%   - mu, Sigma, or alpha: sufficient statistics of the pdf
+%       -> for a Gaussian distribution
+%                mu = E[x]
+%                Sigma = V[x]
+%       -> for a Dirichlet distribution 
+%               alpha = the Dirichlet pseudo-counts
 %
-%       -> for a Dirichlet distribution, set 
-%               moment_1 = the Dirichlet pseudo-counts, alpha
-%          ie, call VBA_exceedanceProbability ('Dirichlet', alpha)
-%
-%   - options: structure with the optional fields:
+%   - options: optional structure with the optional fields:
 %       + verbose: display textual info {false}
 %       + method: 'sampling' or {'analytical'} derivation of the ep
 %       + nSamples: number of samples for the sampling method {5e6}
@@ -80,10 +78,6 @@
     return;
 end
 
-% initialisation
-% =========================================================================
-ep = ones (K, 1);
-
 % ep computation
 % =========================================================================
 
@@ -140,6 +134,10 @@
      
 end
 
+% formating
+% =========================================================================
+ep = VBA_vec (ep);
+
 end
 
 % #########################################################################

From 30d429c3e63e1f8e863b4c8ab4c545e39813c536 Mon Sep 17 00:00:00 2001
From: "lionel.rigoux@sf.mpg.de" <lionel.rigoux@gmail.com>
Date: Mon, 1 Oct 2018 10:11:14 +0200
Subject: [PATCH 14/14] compatibility matlab 2013

---
 utils/VBA_exceedanceProbability.m |  2 +-
 utils/VBA_isdiag.m                | 24 ++++++++++++++++++++++++
 2 files changed, 25 insertions(+), 1 deletion(-)
 create mode 100644 utils/VBA_isdiag.m

diff --git a/utils/VBA_exceedanceProbability.m b/utils/VBA_exceedanceProbability.m
index f9962b47..fd0669ca 100644
--- a/utils/VBA_exceedanceProbability.m
+++ b/utils/VBA_exceedanceProbability.m
@@ -93,7 +93,7 @@
                 ep = ep / sum (ep);
                 
             case 'analytical'
-                if ~ isdiag(Sigma)
+                if ~ VBA_isdiag (Sigma)
                     error('*** VBA_exceedanceProbability: cannot compute analytical ep for non diagonal covariance. Please use the sampling method');
                 end
                 dSigma = max (diag (Sigma), 1e-6);
diff --git a/utils/VBA_isdiag.m b/utils/VBA_isdiag.m
new file mode 100644
index 00000000..e8ffa5b0
--- /dev/null
+++ b/utils/VBA_isdiag.m
@@ -0,0 +1,24 @@
+function flag = VBA_isdiag (A)
+% // VBA toolbox //////////////////////////////////////////////////////////
+%
+% flag = VBA_isdiag (A)
+% check if a matrix is diagonal
+%
+% IN:
+%   - A: matrix
+% OUT:
+%   - flag: boolean flag, true if the matrix A is diagonal
+%
+% Notes:
+% ~~~~~~
+% this is a fallback for isdiag on Matlab < 2014a
+%
+% /////////////////////////////////////////////////////////////////////////
+
+try
+    flag = isdiag (A);
+    return
+catch
+    flag = all (VBA_vec (A == diag (diag (A))));
+end
+