Skip to content

Conversion of Noise Distribution #77

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 11 commits into
base: master
Choose a base branch
from
Open

Conversion of Noise Distribution #77

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
815597a
Add files via upload
May 13, 2019
139985b
Update VBA_Create_NoisePrior.m
May 13, 2019
26c207a
Updating Numerical Optimization
Nov 26, 2019
389f9b2
Error handling
Nov 26, 2019
306dff7
Add files via upload
Nov 27, 2019
222b016
Minor Changes
Nov 29, 2019
e4df934
Delete VBA_Convert_ab.m
Dec 4, 2019
744b3b8
Delete VBA_Create_NoisePrior.m
Dec 4, 2019
60892af
Add files via upload
Dec 4, 2019
b60cd33
Delete demo_dynamicalSystem_Noise.m
Dec 4, 2019
38ca1e2
Add files via upload
Dec 4, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
29 changes: 29 additions & 0 deletions utils/VBA_Convert_ab.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,29 @@
function [mu,sig] = VBA_Convert_ab(a,b)
Copy link
Member

@lionel-rigoux lionel-rigoux Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I find the name a bit pervasive...

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

% [mu,sig] = Convert_ab(a,b)
%
% 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
Copy link
Member

@lionel-rigoux lionel-rigoux Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

%
% M. Eichenlaub 13/05/2019

if isinf(a) % Return 0 if system noise was switched off, i.e. a_sigma = Inf
mu=0;
sig = 0;
elseif a>1
Copy link
Member

@lionel-rigoux lionel-rigoux Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please explain what you are doing here, or provide a reference to an online documentation or notes at the beginning of the function

mu = sqrt(b)*exp(gammaln(a-0.5)-gammaln(a));
sig = sqrt(b/(a-1)-b*exp(2*(gammaln(a-0.5)-gammaln(a))));
else
disp('Conversion error: a<=1');
end
Copy link
Member

@lionel-rigoux lionel-rigoux Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.





50 changes: 50 additions & 0 deletions utils/VBA_Create_NoisePrior.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
function [ a,b ] = VBA_Create_NoisePrior( mu,sig )
Copy link
Member

@lionel-rigoux lionel-rigoux Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

%[ a,b ] = Create_NoisePrior( mu,sig )
%
% 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
Comment on lines +4 to +6
Copy link
Member

@lionel-rigoux lionel-rigoux Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

change to match the doc of the other function


% IN:
% mu,sig: mean and standard deviation of prior distribution over s
Copy link
Member

@lionel-rigoux lionel-rigoux Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

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

@lionel-rigoux lionel-rigoux Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

% M. Eichenlaub 13/05/2019

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

@lionel-rigoux lionel-rigoux Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.


% Starting value
a0 = 1/8*(1+sqrt(49+mu^4/sig^4+50*mu^2/sig^2)+mu^2/sig^2);

% Find a
Copy link
Member

@lionel-rigoux lionel-rigoux Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be nice to have a bit more explanations here about which a you're looking for

options = optimoptions('fmincon','Display','none','Tolfun',1e-8,'TolX',1e-8);
a = fmincon(@(x)D(x,mu,sig),a0,[],[],[],[],1,Inf,[],options);

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

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

@lionel-rigoux lionel-rigoux Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be in the "check results" section


% 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;
Copy link
Member

@lionel-rigoux lionel-rigoux Nov 26, 2019

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)

b=0;
disp('optimization Error');
end

end

end

function [ D ] = D(a,mu,sig)
% Function has to be minimized with respect to a

S = exp(2*(gammaln(a-0.5)-gammaln(a)));
D = mu^2./S - sig^2./((1./(a-1))-S);
D = log(D.^2+1);

end