From 6ac8fb8ffb2e756acf4a54d20fca5e77e6bac3ec Mon Sep 17 00:00:00 2001
From: Daniel Waggoner <dwaggoner@frbatlanta.org>
Date: Fri, 4 May 2012 17:40:33 -0400
Subject: [PATCH] killed bug in dw_elliptical.c from previous commit (cherry
picked from commit 03ce6c37b400f5004a36e0a3104d92bd3d051dbf)
---
elliptical/dw_elliptical.c | 17 ++++++++-------
include/dw_elliptical.h | 42 +-------------------------------------
2 files changed, 9 insertions(+), 50 deletions(-)
diff --git a/elliptical/dw_elliptical.c b/elliptical/dw_elliptical.c
index 9a0501d..db80f84 100644
--- a/elliptical/dw_elliptical.c
+++ b/elliptical/dw_elliptical.c
@@ -153,7 +153,7 @@ static PRECISION draw_truncated_gaussian(PRECISION *draw, TElliptical *elliptica
dw_NormalVector(elliptical->draw);
while ((s=Norm(elliptical->draw)) == 0.0);
r=sqrt(dw_chi_square_invcdf(d->cumulative_lower_bound + dw_uniform_rnd()*d->probability,elliptical->dim));
- ProductVS(elliptical->draw,elliptical->draw,r);
+ ProductVS(elliptical->draw,elliptical->draw,r/s);
ProductMV(elliptical->draw,elliptical->scale,elliptical->draw);
AddVV(elliptical->draw,elliptical->draw,elliptical->center);
if (draw && (draw != pElementV(elliptical->draw)))
@@ -169,19 +169,19 @@ static void print_info_truncated_gaussian(FILE *f_out, TElliptical *elliptical)
}
/*
- Let z be a normal random vector with variance equal to scale*scale' and mean
- equal to center. Truncate z so that it is on the elliptical annulus which is
- given by the set of all z such that
+ Let z be a normal random vector with given center and variance. Truncate z so
+ that it is on the elliptical annulus given by the set of all z such that
- r = sqrt((z - center)' * Inverse(scale * scale') * (z - center))
+ r = sqrt((z - center)' * Inverse(variance) * (z - center))
- is between lower_bound and upper_bound. The density is given by
+ is between lower_bound and upper_bound. Choose scale so that variance=scale*scale'.
+ The density is given by
Gamma(dim/2) f(r)
--------------------------- -----------
2 pi^(dim/2) |det(scale)| r^(dim-1)
- The density of r is given by
+ where the density of r is given by
r^(dim-1)*exp(-(r)^2/2)
@@ -212,8 +212,7 @@ static void print_info_truncated_gaussian(FILE *f_out, TElliptical *elliptical)
The GSL functions gsl_cdf_chisq_P() and gls_cdf_chisq_Pinv() are used to compute
these functions.
*/
-TElliptical* CreateElliptical_TruncatedGaussian(int dim, TVector mean, TMatrix variance,
- PRECISION lower_bound, PRECISION upper_bound)
+TElliptical* CreateElliptical_TruncatedGaussian(int dim, TVector mean, TMatrix variance, PRECISION lower_bound, PRECISION upper_bound)
{
TElliptical* elliptical;
TElliptical_gaussian *data;
diff --git a/include/dw_elliptical.h b/include/dw_elliptical.h
index d70316e..e5a5663 100644
--- a/include/dw_elliptical.h
+++ b/include/dw_elliptical.h
@@ -126,7 +126,7 @@ This gives that the log of the density of z is given by
-ln(abs(det(S)))-ln(2)-0.5*n*ln(pi)+ln(gamma(n/2))-(n-1)*ln(r)+ln(f(r))
-where r=norm(Inverse(S)*(z-c)).
+where r=(z-c)'*Inverse(S*S')*(z-c).
---------------------------------------------------------------------------------
Note that
@@ -153,45 +153,5 @@ The mapping to the data structure is:
quadratic_form = Inverse(S*S')
logabsdet = -ln(abs(det(S)))
dim = n
-
----------------------------------------------------------------------------------
-
-The truncated power case:
-
- Creates a power distribution on the elliptical annulus which is given by the
- set of all z such that
-
- r = sqrt((z - center)' * Inverse(scale * scale') * (z - center))
-
- is between lower_bound and upper_bound. The density is given by
-
- Gamma(dim/2) f(r)
- --------------------------- -----------
- 2 pi^(dim/2) |det(scale)| r^(dim-1)
-
- where f(r) = power * r^(power-1) / (upper_bound^power - lower_bound^power) for
- lower_bound <= r <= upper_bound and power > 0.
-
- Drawing z
- The vector z can be obtained by drawing y from the standard dim-dimensional
- Gaussian distribtuion, r from the distribution on [lower_bound,upper_bound]
- with density f(r), and defining z = r*scale*(y/norm(y)) + center.
-
- Since the cumulative density of r is
-
- (r^power - lower_bound^k) / (upper_bound^power - lower_bound^power)
-
- a draw of r can be obtained by drawing u from the uniform on [0,1] and
- defining by
-
- r = (multiplier*u + shift)^(1/power).
-
- where multiplier = (upper_bound^power - lower_bound^power) and
- shift = lower_bound^power
-
- The variance_scale = sigma^2 * E[r^2] is
- power*(upper_bound^(power+2) - lower_bound^(power+2))
- -------------------------------------------------------
- dim*(power+2)*(upper_bound^power - lower_bound^power)
********************************************************************************/
--
GitLab