From f3c05d6e2caaf3ba975b01efd16c0ab470bd83d8 Mon Sep 17 00:00:00 2001 From: Luigi Acerbi Date: Wed, 26 Oct 2022 18:56:11 +0300 Subject: [PATCH] v1.0.12 Added log pdfs (#6) * feat: added log pdfs and updated example file * fix: updated version * Update README.md --- README.md | 2 +- shared/msmoothboxlogpdf.m | 65 ++++++++++++++ shared/msmoothboxpdf.m | 46 +--------- shared/msplinetrapezlogpdf.m | 71 +++++++++++++++ shared/msplinetrapezpdf.m | 42 +-------- shared/mtrapezlogpdf.m | 63 +++++++++++++ shared/mtrapezpdf.m | 37 +------- shared/munifboxlogpdf.m | 52 +++++++++++ shared/munifboxpdf.m | 27 +----- test/test_pdfs_vbmc.m | 17 ++++ vbmc.m | 6 +- vbmc_examples.m | 165 ++++++++++++++++++++++++----------- 12 files changed, 392 insertions(+), 201 deletions(-) create mode 100644 shared/msmoothboxlogpdf.m create mode 100644 shared/msplinetrapezlogpdf.m create mode 100644 shared/mtrapezlogpdf.m create mode 100644 shared/munifboxlogpdf.m diff --git a/README.md b/README.md index d32ba11..2e21607 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -# Variational Bayesian Monte Carlo (VBMC) - v1.0.11 +# Variational Bayesian Monte Carlo (VBMC) - v1.0.12 **News:** - Added a [Presentations](#presentations) section with links to (relatively) recent slides and video recordings of VBMC related work. diff --git a/shared/msmoothboxlogpdf.m b/shared/msmoothboxlogpdf.m new file mode 100644 index 0000000..3bcb914 --- /dev/null +++ b/shared/msmoothboxlogpdf.m @@ -0,0 +1,65 @@ +function y = msmoothboxlogpdf(x,a,b,sigma) +%MSMOOTHBOXLOGPDF Multivariate smooth-box log probability density function. +% Y = MSMOOTHBOXLOGPDF(X,A,B,SIGMA) returns the logarithm of the pdf of +% the multivariate smooth-box distribution with pivots A and B and scale +% SIGMA, evaluated at the values in X. The multivariate smooth-box pdf is +% the product of univariate smooth-box pdfs in each dimension. +% +% For each dimension i, the univariate smooth-box pdf is defined as a +% uniform distribution between pivots A(i), B(i) and Gaussian tails that +% fall starting from p(A(i)) to the left (resp., p(B(i)) to the right) +% with standard deviation SIGMA(i). +% +% X can be a matrix, where each row is a separate point and each column +% is a different dimension. Similarly, A, B, and SIGMA can also be +% matrices of the same size as X. +% +% The log pdf is typically preferred in numerical computations involving +% probabilities, as it is more stable. +% +% See also MSMOOTHBOXPDF, MSMOOTHBOXRND. + +% Luigi Acerbi 2022 + +[N,D] = size(x); + +if any(sigma(:) <= 0) + error('msmoothboxpdf:NonPositiveSigma', ... + 'All elements of SIGMA should be positive.'); +end + +if D > 1 + if isscalar(a); a = a*ones(1,D); end + if isscalar(b); b = b*ones(1,D); end + if isscalar(sigma); sigma = sigma*ones(1,D); end +end + +if size(a,2) ~= D || size(b,2) ~= D || size(sigma,2) ~= D + error('msmoothboxpdf:SizeError', ... + 'A, B, SIGMA should be scalars or have the same number of columns as X.'); +end + +if size(a,1) == 1; a = repmat(a,[N,1]); end +if size(b,1) == 1; b = repmat(b,[N,1]); end +if size(sigma,1) == 1; sigma = repmat(sigma,[N,1]); end + +if any(a(:) >= b(:)) + error('msmoothboxpdf:OrderError', ... + 'For all elements of A and B, the order A < B should hold.'); +end + +y = -inf(size(x)); +lnf = log(1/sqrt(2*pi)./sigma) - log1p(1/sqrt(2*pi)./sigma.*(b - a)); + +for ii = 1:D + idx = x(:,ii) < a(:,ii); + y(idx,ii) = lnf(idx,ii) - 0.5*((x(idx,ii) - a(idx,ii))./sigma(idx,ii)).^2; + + idx = x(:,ii) >= a(:,ii) & x(:,ii) <= b(:,ii); + y(idx,ii) = lnf(idx,ii); + + idx = x(:,ii) > b(:,ii); + y(idx,ii) = lnf(idx,ii) - 0.5*((x(idx,ii) - b(idx,ii))./sigma(idx,ii)).^2; +end + +y = sum(y,2); \ No newline at end of file diff --git a/shared/msmoothboxpdf.m b/shared/msmoothboxpdf.m index be204ba..558cb56 100644 --- a/shared/msmoothboxpdf.m +++ b/shared/msmoothboxpdf.m @@ -14,50 +14,8 @@ % is a different dimension. Similarly, A, B, and SIGMA can also be % matrices of the same size as X. % -% See also MSMOOTHBOXRND. +% See also MSMOOTHBOXLOGPDF, MSMOOTHBOXRND. % Luigi Acerbi 2022 -[N,D] = size(x); - -if any(sigma(:) <= 0) - error('msmoothboxpdf:NonPositiveSigma', ... - 'All elements of SIGMA should be positive.'); -end - -if D > 1 - if isscalar(a); a = a*ones(1,D); end - if isscalar(b); b = b*ones(1,D); end - if isscalar(sigma); sigma = sigma*ones(1,D); end -end - -if size(a,2) ~= D || size(b,2) ~= D || size(sigma,2) ~= D - error('msmoothboxpdf:SizeError', ... - 'A, B, SIGMA should be scalars or have the same number of columns as X.'); -end - -if size(a,1) == 1; a = repmat(a,[N,1]); end -if size(b,1) == 1; b = repmat(b,[N,1]); end -if size(sigma,1) == 1; sigma = repmat(sigma,[N,1]); end - -if any(a(:) >= b(:)) - error('msmoothboxpdf:OrderError', ... - 'For all elements of A and B, the order A < B should hold.'); -end - -y = zeros(size(x)); -nf = 1 + 1/sqrt(2*pi)./sigma.*(b - a); - -for ii = 1:D - idx = x(:,ii) < a(:,ii); - y(idx,ii) = 1/sqrt(2*pi)./sigma(idx,ii).* exp(-0.5*((x(idx,ii) - a(idx,ii))./sigma(idx,ii)).^2) ./ nf(idx,ii); - - idx = x(:,ii) >= a(:,ii) & x(:,ii) <= b(:,ii); - y(idx,ii) = 1/sqrt(2*pi)./sigma(idx,ii) ./ nf(idx,ii); - - idx = x(:,ii) > b(:,ii); - y(idx,ii) = 1/sqrt(2*pi)./sigma(idx,ii).* exp(-0.5*((x(idx,ii) - b(idx,ii))./sigma(idx,ii)).^2) ./ nf(idx,ii); - -end - -y = prod(y,2); \ No newline at end of file +y = exp(msmoothboxlogpdf(x,a,b,sigma)); \ No newline at end of file diff --git a/shared/msplinetrapezlogpdf.m b/shared/msplinetrapezlogpdf.m new file mode 100644 index 0000000..9bec962 --- /dev/null +++ b/shared/msplinetrapezlogpdf.m @@ -0,0 +1,71 @@ +function y = msplinetrapezlogpdf(x,a,b,c,d) +%MSPLINETRAPEZLOGPDF Multivariate spline-trapezoidal log pdf. +% Y = MSPLINETRAPEZLOGPDF(X,A,B,C,D) returns the logarithm of the pdf of +% the multivariate spline-trapezoidal distribution with external bounds +% A and D and internal points B and C, evaluated at the values in X. The +% multivariate pdf is the product of univariate spline-trapezoidal pdfs +% in each dimension. +% +% For each dimension i, the univariate spline-trapezoidal pdf is defined +% as a trapezoidal pdf whose points A, B and C, D are connected by cubic +% splines such that the pdf is continuous and its derivatives at A, B, C, +% and D are zero (so the derivatives are also continuous): +% +% | __________ +% | /| |\ +% p(X(i)) | / | | \ +% | / | | \ +% |___/___|________|___\____ +% A(i) B(i) C(i) D(i) +% X(i) +% +% X can be a matrix, where each row is a separate point and each column +% is a different dimension. Similarly, A, B, C, and D can also be +% matrices of the same size as X. +% +% The log pdf is typically preferred in numerical computations involving +% probabilities, as it is more stable. +% +% See also MSPLINETRAPEZPDF, MSPLINETRAPEZRND. + +% Luigi Acerbi 2022 + +[N,D] = size(x); + +if D > 1 + if isscalar(a); a = a*ones(1,D); end + if isscalar(b); b = b*ones(1,D); end + if isscalar(c); c = c*ones(1,D); end + if isscalar(d); d = d*ones(1,D); end +end + +if size(a,2) ~= D || size(b,2) ~= D || size(c,2) ~= D || size(d,2) ~= D + error('msplinetrapezlogpdf:SizeError', ... + 'A, B, C, D should be scalars or have the same number of columns as X.'); +end + +if size(a,1) == 1; a = repmat(a,[N,1]); end +if size(b,1) == 1; b = repmat(b,[N,1]); end +if size(c,1) == 1; c = repmat(c,[N,1]); end +if size(d,1) == 1; d = repmat(d,[N,1]); end + +y = -inf(size(x)); +% Normalization factor +% nf = c - b + 0.5*(d - c + b - a); +lnf = log(0.5*(c - b + d - a)); + +for ii = 1:D + idx = x(:,ii) >= a(:,ii) & x(:,ii) < b(:,ii); + z = (x(idx,ii) - a(idx,ii))./(b(idx,ii) - a(idx,ii)); + y(idx,ii) = log(-2*z.^3 + 3*z.^2) - lnf(idx,ii); + + idx = x(:,ii) >= b(:,ii) & x(:,ii) < c(:,ii); + y(idx,ii) = -lnf(idx,ii); + + idx = x(:,ii) >= c(:,ii) & x(:,ii) < d(:,ii); + z = 1 - (x(idx,ii) - c(idx,ii)) ./ (d(idx,ii) - c(idx,ii)); + y(idx,ii) = log(-2*z.^3 + 3*z.^2) - lnf(idx,ii); +end + +y = sum(y,2); + diff --git a/shared/msplinetrapezpdf.m b/shared/msplinetrapezpdf.m index 1a63a15..b984a7e 100644 --- a/shared/msplinetrapezpdf.m +++ b/shared/msplinetrapezpdf.m @@ -22,46 +22,8 @@ % is a different dimension. Similarly, A, B, C, and D can also be % matrices of the same size as X. % -% See also MSPLINETRAPEZRND. +% See also MSPLINETRAPEZLOGPDF, MSPLINETRAPEZRND. % Luigi Acerbi 2022 -[N,D] = size(x); - -if D > 1 - if isscalar(a); a = a*ones(1,D); end - if isscalar(b); b = b*ones(1,D); end - if isscalar(c); c = c*ones(1,D); end - if isscalar(d); d = d*ones(1,D); end -end - -if size(a,2) ~= D || size(b,2) ~= D || size(c,2) ~= D || size(d,2) ~= D - error('msplinetrapezpdf:SizeError', ... - 'A, B, C, D should be scalars or have the same number of columns as X.'); -end - -if size(a,1) == 1; a = repmat(a,[N,1]); end -if size(b,1) == 1; b = repmat(b,[N,1]); end -if size(c,1) == 1; c = repmat(c,[N,1]); end -if size(d,1) == 1; d = repmat(d,[N,1]); end - -y = zeros(size(x)); -% Normalization factor -% nf = c - b + 0.5*(d - c + b - a); -nf = 0.5*(c - b + d - a); - -for ii = 1:D - idx = x(:,ii) >= a(:,ii) & x(:,ii) < b(:,ii); - z = (x(idx,ii) - a(idx,ii))./(b(idx,ii) - a(idx,ii)); - y(idx,ii) = (-2*z.^3 + 3*z.^2) ./ nf(idx,ii); - - idx = x(:,ii) >= b(:,ii) & x(:,ii) < c(:,ii); - y(idx,ii) = 1 ./ nf(idx,ii); - - idx = x(:,ii) >= c(:,ii) & x(:,ii) < d(:,ii); - z = 1 - (x(idx,ii) - c(idx,ii)) ./ (d(idx,ii) - c(idx,ii)); - y(idx,ii) = (-2*z.^3 + 3*z.^2) ./ nf(idx,ii); -end - -y = prod(y,2); - +y = exp(msplinetrapezlogpdf(x,a,b,c,d)); \ No newline at end of file diff --git a/shared/mtrapezlogpdf.m b/shared/mtrapezlogpdf.m new file mode 100644 index 0000000..c6f54bc --- /dev/null +++ b/shared/mtrapezlogpdf.m @@ -0,0 +1,63 @@ +function y = mtrapezlogpdf(x,a,u,v,b) +%MTRAPEZLOGPDF Multivariate trapezoidal probability log pdf. +% Y = MTRAPEZLOGPDF(X,A,U,V,B) returns the logarithm of the pdf of the +% multivariate trapezoidal distribution with external bounds A and B and +% internal points U and V, evaluated at the values in X. The multivariate +% trapezoidal pdf is the product of univariate trapezoidal pdfs in each +% dimension. +% +% For each dimension i, the univariate trapezoidal pdf is defined as: +% +% | __________ +% | /| |\ +% p(X(i)) | / | | \ +% | / | | \ +% |___/___|________|___\____ +% A(i) U(i) V(i) B(i) +% X(i) +% +% X can be a matrix, where each row is a separate point and each column +% is a different dimension. Similarly, A, B, C, and D can also be +% matrices of the same size as X. +% +% The log pdf is typically preferred in numerical computations involving +% probabilities, as it is more stable. +% +% See also MTRAPEZPDF, MTRAPEZRND. + +% Luigi Acerbi 2022 + +[N,D] = size(x); + +if D > 1 + if isscalar(a); a = a*ones(1,D); end + if isscalar(u); u = u*ones(1,D); end + if isscalar(v); v = v*ones(1,D); end + if isscalar(b); b = b*ones(1,D); end +end + +if size(a,2) ~= D || size(u,2) ~= D || size(v,2) ~= D || size(b,2) ~= D + error('mtrapezpdf:SizeError', ... + 'A, B, C, D should be scalars or have the same number of columns as X.'); +end + +if size(a,1) == 1; a = repmat(a,[N,1]); end +if size(u,1) == 1; u = repmat(u,[N,1]); end +if size(v,1) == 1; v = repmat(v,[N,1]); end +if size(b,1) == 1; b = repmat(b,[N,1]); end + +y = -inf(size(x)); +lnf = log(0.5) + log(b - a + v - u) + log(u - a); + +for ii = 1:D + idx = x(:,ii) >= a(:,ii) & x(:,ii) < u(:,ii); + y(idx,ii) = log(x(idx,ii) - a(idx,ii)) - lnf(idx,ii); + + idx = x(:,ii) >= u(:,ii) & x(:,ii) < v(:,ii); + y(idx,ii) = log(u(idx,ii)-a(idx,ii)) - lnf(idx,ii); + + idx = x(:,ii) >= v(:,ii) & x(:,ii) < b(:,ii); + y(idx,ii) = log(b(idx,ii) - x(idx,ii)) - log(b(idx,ii) - v(idx,ii)) + log(u(idx,ii)-a(idx,ii)) - lnf(idx,ii); +end + +y = sum(y,2); \ No newline at end of file diff --git a/shared/mtrapezpdf.m b/shared/mtrapezpdf.m index aa78eff..e006ac7 100644 --- a/shared/mtrapezpdf.m +++ b/shared/mtrapezpdf.m @@ -19,41 +19,8 @@ % is a different dimension. Similarly, A, B, C, and D can also be % matrices of the same size as X. % -% See also MTRAPEZRND. +% See also MTRAPEZLOGPDF, MTRAPEZRND. % Luigi Acerbi 2022 -[N,D] = size(x); - -if D > 1 - if isscalar(a); a = a*ones(1,D); end - if isscalar(u); u = u*ones(1,D); end - if isscalar(v); v = v*ones(1,D); end - if isscalar(b); b = b*ones(1,D); end -end - -if size(a,2) ~= D || size(u,2) ~= D || size(v,2) ~= D || size(b,2) ~= D - error('mtrapezpdf:SizeError', ... - 'A, B, C, D should be scalars or have the same number of columns as X.'); -end - -if size(a,1) == 1; a = repmat(a,[N,1]); end -if size(u,1) == 1; u = repmat(u,[N,1]); end -if size(v,1) == 1; v = repmat(v,[N,1]); end -if size(b,1) == 1; b = repmat(b,[N,1]); end - -y = zeros(size(x)); -nf = 0.5 * (b - a + v - u) .* (u - a); - -for ii = 1:D - idx = x(:,ii) >= a(:,ii) & x(:,ii) < u(:,ii); - y(idx,ii) = (x(idx,ii) - a(idx,ii)) ./ nf(idx,ii); - - idx = x(:,ii) >= u(:,ii) & x(:,ii) < v(:,ii); - y(idx,ii) = (u(idx,ii)-a(idx,ii)) ./ nf(idx,ii); - - idx = x(:,ii) >= v(:,ii) & x(:,ii) < b(:,ii); - y(idx,ii) = (b(idx,ii) - x(idx,ii))./(b(idx,ii) - v(idx,ii)) .* (u(idx,ii)-a(idx,ii)) ./ nf(idx,ii); -end - -y = prod(y,2); \ No newline at end of file +y = exp(mtrapezlogpdf(x,a,u,v,b)); \ No newline at end of file diff --git a/shared/munifboxlogpdf.m b/shared/munifboxlogpdf.m new file mode 100644 index 0000000..5af46a3 --- /dev/null +++ b/shared/munifboxlogpdf.m @@ -0,0 +1,52 @@ +function y = munifboxlogpdf(x,a,b) +%MUNIFBOXLOGPDF Multivariate uniform box log probability density function. +% Y = MUNIFBOXLOGPDF(X,A,B) returns the logarithm of the pdf of the +% multivariate uniform-box distribution with bounds A and B, evaluated at +% the values in X. The multivariate uniform box pdf is the product of +% univariate uniform pdfs in each dimension. +% +% For each dimension i, the univariate uniform-box pdf is defined as: +% +% | +% | ______________ +% p(X(i)) | | | +% | | | +% |___|____________|_____ +% A(i) B(i) +% X(i) +% +% X can be a matrix, where each row is a separate point and each column +% is a different dimension. Similarly, A and B can also be matrices of +% the same size as X. +% +% The log pdf is typically preferred in numerical computations involving +% probabilities, as it is more stable. +% +% See also MUNIFBOXPDF, MUNIFBOXRND. + +% Luigi Acerbi 2022 + +[N,D] = size(x); + +if D > 1 + if isscalar(a); a = a*ones(1,D); end + if isscalar(b); b = b*ones(1,D); end +end + +if size(a,2) ~= D || size(b,2) ~= D + error('munifboxlogpdf:SizeError', ... + 'A, B should be scalars or have the same number of columns as X.'); +end + +if size(a,1) == 1; a = repmat(a,[N,1]); end +if size(b,1) == 1; b = repmat(b,[N,1]); end + +if any(a(:) >= b(:)) + error('munifboxlogpdf:OrderError', ... + 'For all elements of A and B, the order A < B should hold.'); +end + +lnf = sum(log(b - a),2); +y = -lnf .* ones(N,1); +idx = any(bsxfun(@lt, x, a),2) | any(bsxfun(@gt, x, b),2); +y(idx) = -inf; \ No newline at end of file diff --git a/shared/munifboxpdf.m b/shared/munifboxpdf.m index ec2df4b..75ae5fa 100644 --- a/shared/munifboxpdf.m +++ b/shared/munifboxpdf.m @@ -19,31 +19,8 @@ % is a different dimension. Similarly, A and B can also be matrices of % the same size as X. % -% See also MUNIFBOXRND. +% See also MUNIFBOXLOGPDF, MUNIFBOXRND. % Luigi Acerbi 2022 -[N,D] = size(x); - -if D > 1 - if isscalar(a); a = a*ones(1,D); end - if isscalar(b); b = b*ones(1,D); end -end - -if size(a,2) ~= D || size(b,2) ~= D - error('munifboxpdf:SizeError', ... - 'A, B should be scalars or have the same number of columns as X.'); -end - -if size(a,1) == 1; a = repmat(a,[N,1]); end -if size(b,1) == 1; b = repmat(b,[N,1]); end - -if any(a(:) >= b(:)) - error('munifboxpdf:OrderError', ... - 'For all elements of A and B, the order A < B should hold.'); -end - -nf = prod(b - a,2); -y = 1 ./ nf .* ones(N,1); -idx = any(bsxfun(@lt, x, a),2) | any(bsxfun(@gt, x, b),2); -y(idx) = 0; \ No newline at end of file +y = exp(munifboxlogpdf(x,a,b)); \ No newline at end of file diff --git a/test/test_pdfs_vbmc.m b/test/test_pdfs_vbmc.m index 330ae54..94c3f4d 100644 --- a/test/test_pdfs_vbmc.m +++ b/test/test_pdfs_vbmc.m @@ -13,10 +13,15 @@ function test_pdfs_vbmc() %% Test multivariate uniform box distribution pdf1 = @(x) munifboxpdf(x,a(1),b(1)); pdf2 = @(x) munifboxpdf(x,a,b); +pdf1log = @(x) exp(munifboxlogpdf(x,a(1),b(1))); +pdf2log = @(x) exp(munifboxlogpdf(x,a,b)); pdfrnd = @(n) munifboxrnd(a,b,n); name = 'munifbox'; + test_pdf1_normalization(pdf1,lb,ub,tolerr,name); test_pdf2_normalization(pdf2,lb,ub,tolerr,name); +test_pdf1_normalization(pdf1log,lb,ub,tolerr,name); +test_pdf2_normalization(pdf2log,lb,ub,tolerr,name); test_rnd(pdfrnd,pdf1,a,b,n,tolrmse,name); %% Test multivariate trapezoidal distribution @@ -24,31 +29,43 @@ function test_pdfs_vbmc() v = [1.5,-3.4]; pdf1 = @(x) mtrapezpdf(x,a(1),u(1),v(1),b(1)); pdf2 = @(x) mtrapezpdf(x,a,u,v,b); +pdf1log = @(x) exp(mtrapezlogpdf(x,a(1),u(1),v(1),b(1))); +pdf2log = @(x) exp(mtrapezlogpdf(x,a,u,v,b)); pdfrnd = @(n) mtrapezrnd(a,u,v,b,n); name = 'mtrapez'; test_pdf1_normalization(pdf1,lb,ub,tolerr,name); test_pdf2_normalization(pdf2,lb,ub,tolerr,name); +test_pdf1_normalization(pdf1log,lb,ub,tolerr,name); +test_pdf2_normalization(pdf2log,lb,ub,tolerr,name); test_rnd(pdfrnd,pdf1,a,b,n,tolrmse,name); %% Test multivariate spline trapezoidal distribution pdf1 = @(x) msplinetrapezpdf(x,a(1),u(1),v(1),b(1)); pdf2 = @(x) msplinetrapezpdf(x,a,u,v,b); +pdf1log = @(x) exp(msplinetrapezlogpdf(x,a(1),u(1),v(1),b(1))); +pdf2log = @(x) exp(msplinetrapezlogpdf(x,a,u,v,b)); pdfrnd = @(n) msplinetrapezrnd(a,u,v,b,n); name = 'msplinetrapez'; test_pdf1_normalization(pdf1,lb,ub,tolerr,name); test_pdf2_normalization(pdf2,lb,ub,tolerr,name); +test_pdf1_normalization(pdf1log,lb,ub,tolerr,name); +test_pdf2_normalization(pdf2log,lb,ub,tolerr,name); test_rnd(pdfrnd,pdf1,a,b,n,tolrmse,name); %% Test multivariate smoothbox distribution sigma = [0.7,0.45]; pdf1 = @(x) msmoothboxpdf(x,a(1),b(1),sigma(1)); pdf2 = @(x) msmoothboxpdf(x,a,b,sigma); +pdf1log = @(x) exp(msmoothboxlogpdf(x,a(1),b(1),sigma(1))); +pdf2log = @(x) exp(msmoothboxlogpdf(x,a,b,sigma)); pdfrnd = @(n) msmoothboxrnd(a,b,sigma,n); name = 'msmoothbox'; lb = [-5,-7]; ub = [5,0]; test_pdf1_normalization(pdf1,lb,ub,tolerr,name); test_pdf2_normalization(pdf2,lb,ub,tolerr,name); +test_pdf1_normalization(pdf1log,lb,ub,tolerr,name); +test_pdf2_normalization(pdf2log,lb,ub,tolerr,name); test_rnd(pdfrnd,pdf1,a,b,n,tolrmse,name); %close all; diff --git a/vbmc.m b/vbmc.m index 705df6b..99d1d1e 100644 --- a/vbmc.m +++ b/vbmc.m @@ -1,6 +1,6 @@ function [vp,elbo,elbo_sd,exitflag,output,samples,optimState,stats,vp_train] = ... vbmc(fun,x0,LB,UB,PLB,PUB,options,varargin) -%VBMC Posterior and model inference via Variational Bayesian Monte Carlo (v1.0.11). +%VBMC Posterior and model inference via Variational Bayesian Monte Carlo (v1.0.12). % VBMC computes a variational approximation of the full posterior and a % lower bound on the normalization constant (marginal likelhood or model % evidence) for a provided unnormalized log posterior. As of v1.0, VBMC @@ -140,7 +140,7 @@ % Author (copyright): Luigi Acerbi, 2018-2022 % e-mail: luigi.acerbi@{helsinki.fi,gmail.com} % URL: http://luigiacerbi.com -% Version: 1.0.11 +% Version: 1.0.12 % Release date: Oct 26, 2022 % Code repository: https://github.com/lacerbi/vbmc %-------------------------------------------------------------------------- @@ -149,7 +149,7 @@ %% Start timer t0 = tic; -vbmc_version = '1.0.11'; +vbmc_version = '1.0.12'; % Check that all VBMC subfolders are on the MATLAB path add2path(); diff --git a/vbmc_examples.m b/vbmc_examples.m index 091eab3..ec6a3f3 100644 --- a/vbmc_examples.m +++ b/vbmc_examples.m @@ -472,50 +472,66 @@ fprintf(' *** Bounded parameters:\n\n'); fprintf(' The first simple choice is a non-informative uniformly flat prior between the hard bounds.\n'); fprintf(' This is implemented in the `munifboxpdf` (multivariate uniform-box) function.\n'); -fprintf(' We plot the marginal pdf of the prior for two-dimensional problem (see plot, 1st row).\n\n'); +fprintf(' We plot the pdf of the prior for a one-dimensional problem (see plot, 1st row),\n'); +fprintf(' the prior pdf to the left and the log pdf, as requested by VBMC, to the right\n'); +fprintf(' (implemented by the `munifboxlogpdf` function).\n\n'); fprintf(' Press any key to continue.\n\n'); -D = 2; -lb = [-3,-6]; -ub = [3,6]; -plb = [-2,-3]; -pub = [2,3]; +lb = -3; +ub = 3; +plb = -2; +pub = 2; -parameter_names{1} = 'x_1'; -parameter_names{2} = 'x_2'; -ylims = [0,0.25]; +ylims = [0,0.25; -20, 0]; figure(1); -for d = 1:D - x = linspace(lb(d)-1,ub(d)+1,1e3)'; - subplot(3,D,d); - plot(x, munifboxpdf(x,lb(d),ub(d)), 'k-', 'LineWidth', 1); +for i = 1:2 + x = linspace(lb-1,ub+1,1e3)'; + switch i + case 1; y = munifboxpdf(x,lb,ub); + case 2; y = munifboxlogpdf(x,lb,ub); + end + subplot(3,2,i); + plot(x, max(y,ylims(i,1)), 'k-', 'LineWidth', 1); set(gca,'TickDir','out'); - xlabel(parameter_names{d}) - ylabel('prior pdf') - ylim(ylims); + xlabel('x_1') + switch i + case 1; ylabel('prior pdf'); + case 2; ylabel('prior log pdf'); + end + xlim([x(1),x(end)]); + ylim(ylims(i,:)); box off; - if d == 1; title('Uniform prior'); end + if i == 1; title('Uniform prior'); end end set(gcf,'Color','w'); pause; fprintf(' Alternatively, for each parameter we can define the prior to be flat within a range,\n'); fprintf(' where a reasonable choice is the "plausible" range, and then falls to zero towards the hard bounds.\n'); -fprintf(' This is a trapezoidal or "tent" prior, implemented by the provided `mtrapezpdf` function (see plot, 2nd row).\n\n'); +fprintf(' This is a trapezoidal or "tent" prior, implemented by the provided `mtrapezpdf` and `mtrapezlogpdf`\n'); +fprintf(' functions (see plot, 2nd row).\n\n'); fprintf(' Press any key to continue.\n\n'); figure(1); -for d = 1:D - x = linspace(lb(d)-1,ub(d)+1,1e3)'; - subplot(3,D,d+D); - plot(x, mtrapezpdf(x,lb(d),plb(d),pub(d),ub(d)), 'k-', 'LineWidth', 1); +for i = 1:2 + x = linspace(lb-1,ub+1,1e3)'; + switch i + case 1; y = mtrapezpdf(x,lb,plb,pub,ub); + case 2; y = mtrapezlogpdf(x,lb,plb,pub,ub); + end + subplot(3,2,i+2); + plot(x, max(y,ylims(i,1)), 'k-', 'LineWidth', 1); set(gca,'TickDir','out'); - xlabel(parameter_names{d}) - ylabel('prior pdf') - ylim(ylims); + xlabel('x_1') + switch i + case 1; ylabel('prior pdf'); + case 2; ylabel('prior log pdf'); + end + xlim([x(1),x(end)]); + ylim(ylims(i,:)); box off; - if d == 1; title('Trapezoidal prior'); end + if i == 1; title('Trapezoidal prior'); end end set(gcf,'Color','w'); pause; @@ -524,20 +540,29 @@ fprintf(' (obtained using cubic splines). This prior is better behaved numerically\n'); fprintf(' as it is continuous with continuous derivatives (i.e., no sharp edges),\n'); fprintf(' so we recommend it over the simple trapezoidal prior.\n'); -fprintf(' The spline-smoothed trapezoidal prior is implemented in the `msplinetrapezpdf` function (see plot, 3rd row).\n\n'); +fprintf(' The spline-smoothed trapezoidal prior is implemented in the `msplinetrapezpdf` and\n'); +fprintf(' `msplinetrapezlogpdf` functions (see plot, 3rd row).\n\n'); fprintf(' Press any key to continue.\n\n'); figure(1); -for d = 1:D - x = linspace(lb(d)-1,ub(d)+1,1e3)'; - subplot(3,D,d+2*D); - plot(x, msplinetrapezpdf(x,lb(d),plb(d),pub(d),ub(d)), 'k-', 'LineWidth', 1); +for i = 1:2 + x = linspace(lb-1,ub+1,1e3)'; + switch i + case 1; y = msplinetrapezpdf(x,lb,plb,pub,ub); + case 2; y = msplinetrapezlogpdf(x,lb,plb,pub,ub); + end + subplot(3,2,i+4); + plot(x, max(y,ylims(i,1)), 'k-', 'LineWidth', 1); set(gca,'TickDir','out'); - xlabel(parameter_names{d}) - ylabel('prior pdf') - ylim(ylims); + xlabel('x_1') + switch i + case 1; ylabel('prior pdf'); + case 2; ylabel('prior log pdf'); + end + xlim([x(1),x(end)]); + ylim(ylims(i,:)); box off; - if d == 1; title('Smoothed trapezoidal prior'); end + if i == 1; title('Smoothed trapezoidal prior'); end end set(gcf,'Color','w'); pause; @@ -549,13 +574,13 @@ fprintf(' As an alternative to these common choices, we also provide a smoothbox distribution\n'); fprintf(' which is uniform within an interval (typically, the plausible range) and then falls with\n'); fprintf(' Gaussian tails with scale sigma outside the interval.\n'); +fprintf(' Find the implementation in `msmoothboxpdf` and `msmoothboxlogpdf`.\n'); close all; -D = 2; -lb = -inf(1,D); -ub = inf(1,D); -plb = [0,-2]; -pub = [5,1]; +lb = -inf; +ub = inf; +plb = -1; +pub = 3; % We recommend setting sigma as a fraction of the plausible range. % For example sigma set to 4/10 of the plausible range assigns ~50% @@ -568,23 +593,55 @@ sigma = 0.4*prange; figure(1); -for d = 1:D - subplot(1,D,d); - x = linspace(plb(d) - 2*prange(d), pub(d) + 2*prange(d), 1e3)'; - y = msmoothboxpdf(x, plb(d), pub(d), sigma(d)); - plot(x, y, 'k-', 'LineWidth', 1); +for i = 1:2 + x = linspace(plb - 2*prange, pub + 2*prange, 1e3)'; + switch i + case 1; y = msmoothboxpdf(x, plb, pub, sigma); + case 2; y = msmoothboxlogpdf(x, plb, pub, sigma); + end + subplot(1,2,i); + plot(x, max(y,ylims(i,1)), 'k-', 'LineWidth', 1); set(gca,'TickDir','out'); - xlabel(parameter_names{d}) - ylabel('prior pdf') - ylim(ylims); + xlabel('x_1') + switch i + case 1; ylabel('prior pdf'); + case 2; ylabel('prior log pdf'); + end + xlim([x(1),x(end)]); + ylim(ylims(i,:)); box off; - if d == 1; title('Smoothed box prior (unbounded parameters)'); end + if i == 1; title('Smoothed box prior (unbounded parameters)'); end end set(gcf,'Color','w'); fprintf('\n Press any key to continue.\n\n'); pause; +fprintf(' To conclude, we rerun inference on a bounded problem, using one of the priors\n'); +fprintf(' (smoothed trapezoidal) introduced in this section. See code for details.\n\n'); +fprintf(' Press any key to continue.\n\n'); +pause; + +D = 2; % Still in 2-D +LB = zeros(1,D); +UB = 10*ones(1,D); +PLB = 0.1*ones(1,D); +PUB = 3*ones(1,D); + +% Define the log prior and log likelihood +lpriorfun = @(x) msplinetrapezlogpdf(x,LB,PLB,PUB,UB); +llfun = @rosenbrock_test; +fun = @(x) llfun(x) + lpriorfun(x); + +x0 = ones(1,D); +options = vbmc('defaults'); +[vp,elbo,elbo_sd] = vbmc(fun,x0,LB,UB,PLB,PUB,options); +Xs = vbmc_rnd(vp,3e5); +cornerplot(Xs,{'x_1','x_2'}); + +fprintf(' Press any key to continue.\n\n'); +pause; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Example 6: Noisy log-likelihood evaluations @@ -639,10 +696,12 @@ D = 3; -% We use a smoothed "tent" or trapezoidal prior over the parameters. This -% prior is flat between PLB and PUB, and decreases smoothly to 0 towards LB -% and UB. The prior is implemented using cubic splines, hence the name. -lpriorfun = @(x) log(msplinetrapezpdf(x,LB,PLB,PUB,UB)); +% We use a smoothed "tent" or trapezoidal prior over the parameters (see +% Example 5 in this file). This prior is flat between PLB and PUB, and +% decreases smoothly to 0 towards LB and UB. The prior is implemented +% using cubic splines, hence the name. Note that we pass the function that +% directly computes the log pdf. +lpriorfun = @(x) msplinetrapezlogpdf(x,LB,PLB,PUB,UB); % Since LLFUN has now two outputs (the log likelihood at X, and the % estimated SD of the log likelihood at X), we cannot directly sum the log