blob: bfde566386cd60292679e58e20dc806c7cbd1360 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
|
(sc-include "macros/macros")
(pre-include "gsl/gsl_blas.h")
(pre-include "gsl/gsl_integration.h")
(pre-include "gsl/gsl_randist.h")
(pre-include "utils.h")
(pre-define INTEGRATION-EPSABS 1e-3)
(pre-define INTEGRATION-EPSREL 1e-3)
(define (shifted-gaussian-random-vector r mean theta-max standard-deviation x)
((unsigned int) (const gsl-rng*) (const gsl-vector*) double double gsl-vector*)
"Return a random vector generated by the shifted Gaussian method."
(for-i i (: x size)
(gsl-vector-set x i (gsl-ran-gaussian r standard-deviation)))
(gsl-vector-add x mean)
(gsl-blas-dscal (/ 1 (gsl-blas-dnrm2 x)) x)
(cond
((> (angle-between-vectors x mean) theta-max)
(return (+ 1 (shifted-gaussian-random-vector r mean theta-max standard-deviation x))))
(else (return 1))))
(define (shifted-gaussian-pdf theta mean theta-max standard-deviation dimension w)
(double double double double double (unsigned int) gsl-integration-workspace*)
"Return the unnormalized probability density function at THETA for a
random vector generated using the shifted Gaussian method."
(cond
((> theta theta-max) (return 0))
(else
(define projection double (* mean (cos theta)))
(define (integrand r params) (double double void*)
(return (* (gsl-pow-uint r (- dimension 1))
(gaussian-pdf (/ (- r projection)
standard-deviation)))))
(declare gsl-integrand (struct-variable gsl-function (function (address-of integrand)))
result double
error double)
;; (gsl-integration-qagiu (address-of gsl-integrand)
;; 0
;; INTEGRATION-EPSABS
;; INTEGRATION-EPSREL
;; (: w limit)
;; w
;; (address-of result)
;; (address-of error))
(let* ((rmax double (/ (+ projection
(sqrt (+ (gsl-pow-2 projection)
(* 4
(- dimension 1)
(gsl-pow-2 standard-deviation)))))
2))
(half-width double (* 3 standard-deviation)))
(gsl-integration-qag (address-of gsl-integrand)
(if* (< (- rmax half-width) 0)
0
(- rmax half-width))
(+ rmax half-width)
INTEGRATION-EPSABS
INTEGRATION-EPSREL
(: w limit)
GSL-INTEG-GAUSS41
w
(address-of result)
(address-of error)))
(return (/ (* result
(gaussian-pdf (/ (sqrt (- (gsl-pow-2 mean)
(gsl-pow-2 projection)))
standard-deviation)))
(gsl-pow-uint (sqrt (* 2 M-PI)) (- dimension 2))
(gsl-pow-uint standard-deviation dimension))))))
|
