-
Notifications
You must be signed in to change notification settings - Fork 40
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
* feat: added log pdfs and updated example file * fix: updated version * Update README.md
- Loading branch information
Showing
12 changed files
with
392 additions
and
201 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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); |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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); | ||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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); |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.