Conversion of Noise Distribution #77

manueich · 2019-05-13T13:50:27Z

The two functions added allow the conversion between a gamma distributed precision and the distribution over the associated standard deviation in both directions for the more intuitive interpretation/definition of posterior/prior distributions of the noise parameters.
The attached document provides the mathematical details.
Notes.pdf

lionel-rigoux · 2019-05-13T14:21:07Z

Thank you for your contribution! This looks like very promising and a welcome addition to the toolbox (the current VBA_guess_hyperpriors is a bit rough...). Jean and I will review your code and notes and we'll come back to you if we have questions.
Also, in the case we would integrate your code, would you be willing to write a litte paragraph/page for our wiki to explain how to use your function (for the less experienced users)?

manueich · 2019-05-13T14:39:19Z

Thanks for your feedback. I put it up after Jean suggested it in recent forum post I made. I would be happy to write something for the wiki, should it be accepted.
Something I should mention is that the conversion from mu and sig to a and b is based on a numerical optimization scheme, which does not work for all possible combinations of mu and sig. For example if mu is very large in comparison to sig, e.g. mu=10000 and sig=0.01, the optimization fails. However these cases should be quite unrealistic in practice.
The function automatically checks if the optimization worked within a certain error margin and alerts the user if it failed.

Looking forward to hearing from you.

manueich · 2019-08-23T15:22:21Z

Dear Jean and Lionel,
I was wondering whether you had a chance to review my submission. I have been developing this method in the context of my PhD project and it has fed into a few papers (one of them is listed in VBA's website list of publications under Eichenlaub M. and more are likely to follow), so I would be really interested in your opinion on this. It is not a big deal, but I said I think this is useful to anyone who wants a more intuitive interpretation of the noise parameter distributions.
Best wishes,
Manuel

manueich · 2019-11-24T19:07:33Z

Hi,
I was wondering if there is any chance that this will get reviewed at some point.
Manuel

lionel-rigoux · 2019-11-25T15:13:26Z

Hi Manuel,

Sorry for the long silence, it's been a busy year and the toolbox was quite down the priority list unfortunately.

I'll add some comments to your pull request right away.

The only thing I think is missing is a demo (like demo_static or demo_multisession) with a very simple model that shows the difference between using the default method and the one you suggest. The goal is double: 1) provide the other users with a minimal working example that could serve as a tutorial, 2) have a formal way to check that your method indeed works (either by running it manually, but also through VBA_test). Would you have time to implement this?

Could you also send me a little text for the wiki?

lionel-rigoux

My main points are:

need more doc about what this exactly does (approximation) and what are the limit cases when this function should not be used
proper argument checks and error throwing

I would also love to check the demo, I am curious to see how much better is the parameter recovery using the neat trick. I would also like to compare your approach to my no-so-neat VBA_guessHyperpriors.

lionel-rigoux · 2019-11-26T08:14:06Z

utils/VBA_Convert_ab.m

@@ -0,0 +1,29 @@
+function [mu,sig] = VBA_Convert_ab(a,b)


I find the name a bit pervasive...

Suggested change

function [mu,sig] = VBA_Convert_ab(a,b)

% function [mu, sig] = VBA_Gam2Norm (a, b)

lionel-rigoux · 2019-11-26T08:15:42Z

utils/VBA_Create_NoisePrior.m

@@ -0,0 +1,50 @@
+function [ a,b ] = VBA_Create_NoisePrior( mu,sig )


Suggested change

function [ a,b ] = VBA_Create_NoisePrior( mu,sig )

function [a, b] = VBA_Norm2Gam (mu , sig)```

lionel-rigoux · 2019-11-26T08:19:12Z

utils/VBA_Create_NoisePrior.m

+if nargin==1 % Return this if system noise shall be switched off
+    a = Inf;
+    b = 0;
+else


I would rather enforce having two parameters, return [Inf 0] if sigma is 0.
What happen if the user provides a negative mu or sig? Proper checks need to be done here, with meaningful error messages if they don't pass.

lionel-rigoux · 2019-11-26T08:19:49Z

utils/VBA_Create_NoisePrior.m

+% defined by a and b with p = 1/s^2
+
+% IN:
+%   mu,sig: mean and standard deviation of prior distribution over s


Normal priors are usually defined in the toolbox by their mean and variance (sigma2). It shouldn't be difficult to change your code so it doesn't confuse too much the users.

lionel-rigoux · 2019-11-26T08:24:43Z

utils/VBA_Convert_ab.m

+else
+    disp('Conversion error: a<=1');
+end


This should be a separate if at the beginning which returns an error (with a message starting with *** VBA_Gam2Norm: ) providing details about why this is wrong and how the user can mitigate the issue.

lionel-rigoux · 2019-11-26T08:29:44Z

utils/VBA_Create_NoisePrior.m

+
+b = (mu./(exp(gammaln(a-0.5)-gammaln(a)))).^2;
+
+[m,s] = VBA_Convert_ab(a,b);


This should be in the "check results" section

lionel-rigoux · 2019-11-26T08:40:09Z

utils/VBA_Create_NoisePrior.m

+
+% Check results, reject if error on either mu or sig is greater than 1%
+if abs(m-mu)/mu>1e-2 || abs(s-sig)/mu>1e-2
+    a=0;


This is dangerous. If the approximation fails, this message could be lost in the logs and it would be difficult to detect where the issue comes from (or even to detect that those are weird values for a lot of users).
This case should return an error (see my other comment)

lionel-rigoux · 2019-11-26T08:52:06Z

utils/VBA_Convert_ab.m

+% Converts a and b defining a gamma distributed precision p to mean (mu) and
+% standard deviation (sig) of the distribution over the associated standard
+% deviation s, i.e. s = 1/sqrt(p)
+% IN:
+%   a,b: Can either be a_sigma/b_sigma or a_alpha/b_alpha defining
+%   posterior distributions over measurement/system noise precision
+% OUT:
+%   mu,sig: mean and standard deviation of distribution over s
+% 
+% Be aware that this is only possible for a>1


Convert a Gamma distribution -- defined by the shape and rate parameters a and b -- over the noise precision s into a Normal distribution -- defined by the mean and variance mu and sigma2 -- over the noise standard deviation (1/sqrt(s)).
IN:
a, b: shape and rate of the Gamma distribution over the noise precision (eg. a_sigma and b_sigma)
OUT:
mu, sigma2: mean and variance of the corresponding normal distribution over the noise standard deviation

Note that this transformation is only possible for a > 1.
% Please explain here that this is not an exact transformation but an approximate mapping obtained my matching moments. Provide limit cases when this mapping is likely to fail.

lionel-rigoux · 2019-11-26T08:58:28Z

utils/VBA_Create_NoisePrior.m

+% Converts mu and sig defining a prior distribution over noise standardard
+% deviation s, to a gamma distriribution over the associated precision p
+% defined by a and b with p = 1/s^2


change to match the doc of the other function

lionel-rigoux · 2019-11-26T08:59:09Z

utils/VBA_Create_NoisePrior.m

+% OUT:
+%   a,b: Can either be a_sigma/b_sigma or a_alpha/b_alpha defining
+%   prior distributions over measurement/system noise precision
+% 


as above, it should be explicit that this is only a matching moment approximation of the problem

manueich · 2019-11-26T15:39:53Z

Hello Lionel,
thank you very much for your reply and comments.
Before we go any further, I think I need to clarify something.

Reading your comments, I think you might have slightly misunderstood what the functions are implementing. It seems to me that you think the function does the following: Given a gamma distributed precision p, defined by a and b, the function VBA_Convert_ab (or VBA_Gam2Norm as you suggested) approximates the distribution over the standard deviation s=1/sqrt(p) with a normal distribution, defined with mu and sigma.

This is not exactly the case as I don't use a normal distribution at all. What I suggest is not an approximation, it is the exact thing. Given the transformation s=1/sqrt(p), it is possible to directly formulate the PDF over s, called IRGamma defined by a and b from the "original" gamma distribution. Subsequently, mean and standard deviation, which I call mu and sig, of this new PDF (IRGamma) can be directly calculated using a and b. This is what the function VBA_Convert_ab does.

The second function (VBA_Create_Noise _Prior) essentially reverses this by calculating a and b given mu and sig. Here I use a numerical approximation simply because a direct analytical derivation of a and b is not possible (as far as I can tell and I have tried a few things).

The document I attached in my very first post, provides the details of this. Nevertheless I attached a very short document explaining this as well.

Short.pdf

If my reasoning is correct and you still think this is a valuable contribution to the toolbox, I am happy to make the changes you asked for and write a demo and something for the wiki. I am also happy to test the function against the existing functions if you give me a bit more details on which to use. I already had a look at the VBA_guessHyperpriors function, but am not quite sure how it works, maybe you could give me an example.

Manuel

VBA_Create_NoisePrior: I have now established bounds on mu and sigma based on numerical simulations. These are chosen conservatively and should include all reasonable values occurring in practice. - mu > 5e-3 - 0.005*mu < sig < 5*mu Respective error messages are returned if mu/sigma lie outside these bounds.

Demo on how to use the functions and comparison to VBA_guessHyperpriors

lionel-rigoux · 2019-12-04T10:10:14Z

Hi Manu,

Of course you're not relying on a Gaussian approximation, I don't know why I was confused.

You demo shows well how your function can help the estimation. I'm also glad to see that mine was not too bad either (compared to the default).

I think that I will go one step further and create a more generic function allowing to define the hyperpriors using different summary statistics, like mean and variance of the std (using your function), confidence interval about the explained variance (like in guessHyperpriors), and so on.
To do so, I will create a new branch in which I will accept your pull request. Once the feature is done it will be merged in the dev then master branch. I hope this is fine with you. If you want to help, you could eg. see if deriving a and b from a confidence interval about the noise std is possible.

I just have one request left: you pushed your functions twice, once in the root folder, and once in the utils folder. I am not sure those two are the same... Please delete the outdated one so I am sure to integrate the right one.

New code is working and will be refactored in a feat branch anyway

manueich · 2019-12-04T12:52:12Z

Hi Lionel,

thank you very much. I am happy to hear that you think my contributions are useful. I think the idea of creating a more generic function combining different approaches is great and I look forward to using it.

I apologise for the confusion with the uploads. The two relevant and up to date functions are now only in the utils folder and I have moved the demo into demos/0_basics folder.

Please let me know if I can help with anything else.

Best wishes,
Manuel

manueich · 2021-02-02T13:36:13Z

Hi,

just to give an update on this. I published a preprint on this method (https://arxiv.org/abs/2101.06289) and you will find an updated version of the code in the repository (https://github.com/manueich/Noise_Gamma), in case you decide to include it.

Manuel

Add files via upload

815597a

Update VBA_Create_NoisePrior.m

139985b

lionel-rigoux previously requested changes Nov 26, 2019

View reviewed changes

manueich added 4 commits November 26, 2019 16:03

Updating Numerical Optimization

26c207a

Add files via upload

306dff7

Demo on how to use the functions and comparison to VBA_guessHyperpriors

Minor Changes

222b016

manueich added 5 commits December 4, 2019 12:43

Delete VBA_Convert_ab.m

e4df934

Delete VBA_Create_NoisePrior.m

744b3b8

Add files via upload

60892af

Delete demo_dynamicalSystem_Noise.m

b60cd33

Add files via upload

38ca1e2

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Conversion of Noise Distribution #77

Conversion of Noise Distribution #77

manueich commented May 13, 2019

lionel-rigoux commented May 13, 2019

manueich commented May 13, 2019

manueich commented Aug 23, 2019

manueich commented Nov 24, 2019

lionel-rigoux commented Nov 25, 2019 •

edited

Loading

lionel-rigoux left a comment

lionel-rigoux Nov 26, 2019

lionel-rigoux Nov 26, 2019

lionel-rigoux Nov 26, 2019

lionel-rigoux Nov 26, 2019

lionel-rigoux Nov 26, 2019

lionel-rigoux Nov 26, 2019

lionel-rigoux Nov 26, 2019

lionel-rigoux Nov 26, 2019

lionel-rigoux Nov 26, 2019

lionel-rigoux Nov 26, 2019

manueich commented Nov 26, 2019

lionel-rigoux commented Dec 4, 2019

manueich commented Dec 4, 2019

manueich commented Feb 2, 2021

	function [mu,sig] = VBA_Convert_ab(a,b)
	% function [mu, sig] = VBA_Gam2Norm (a, b)

		@@ -0,0 +1,50 @@
		function [ a,b ] = VBA_Create_NoisePrior( mu,sig )

	function [ a,b ] = VBA_Create_NoisePrior( mu,sig )
	function [a, b] = VBA_Norm2Gam (mu , sig)```


		b = (mu./(exp(gammaln(a-0.5)-gammaln(a)))).^2;

		[m,s] = VBA_Convert_ab(a,b);

Conversion of Noise Distribution #77

Are you sure you want to change the base?

Conversion of Noise Distribution #77

Conversation

manueich commented May 13, 2019

lionel-rigoux commented May 13, 2019

manueich commented May 13, 2019

manueich commented Aug 23, 2019

manueich commented Nov 24, 2019

lionel-rigoux commented Nov 25, 2019 • edited Loading

lionel-rigoux left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

manueich commented Nov 26, 2019

lionel-rigoux commented Dec 4, 2019

manueich commented Dec 4, 2019

manueich commented Feb 2, 2021

lionel-rigoux commented Nov 25, 2019 •

edited

Loading