Commit fb2c19e
| Original file line number | Diff line number | Diff line change | |
|---|---|---|---|
@@ -0,0 +1,133 @@ | |||
| 1 | + module random | ||
| 2 | + | ||
| 3 | + use types, only: dp | ||
| 4 | + use utils, only: stop_error | ||
| 5 | + implicit none | ||
| 6 | + private | ||
| 7 | + public randn, rand_gamma | ||
| 8 | + | ||
| 9 | + interface randn | ||
| 10 | + module procedure randn_scalar | ||
| 11 | + module procedure randn_vector | ||
| 12 | + module procedure randn_matrix | ||
| 13 | + module procedure randn_vector_n | ||
| 14 | + end interface | ||
| 15 | + | ||
| 16 | + interface rand_gamma | ||
| 17 | + module procedure rand_gamma_scalar | ||
| 18 | + module procedure rand_gamma_vector | ||
| 19 | + module procedure rand_gamma_matrix | ||
| 20 | + module procedure rand_gamma_vector_n | ||
| 21 | + end interface | ||
| 22 | + | ||
| 23 | + contains | ||
| 24 | + | ||
| 25 | + subroutine randn_scalar(x) | ||
| 26 | + ! Returns a psuedorandom scalar drawn from the standard normal distribution. | ||
| 27 | + ! | ||
| 28 | + ! [1] Marsaglia, G., & Bray, T. A. (1964). A Convenient Method for Generating | ||
| 29 | + ! Normal Variables. SIAM Review, 6(3), 260–264. | ||
| 30 | + real(dp), intent(out) :: x | ||
| 31 | + logical, save :: first = .true. | ||
| 32 | + real(dp), save :: u(2) | ||
| 33 | + real(dp) :: r2 | ||
| 34 | + if (first) then | ||
| 35 | + do | ||
| 36 | + call random_number(u) | ||
| 37 | + u = 2*u-1 | ||
| 38 | + r2 = sum(u**2) | ||
| 39 | + if (r2 < 1 .and. r2 > 0) exit | ||
| 40 | + end do | ||
| 41 | + u = u * sqrt(-2*log(r2)/r2) | ||
| 42 | + x = u(1) | ||
| 43 | + else | ||
| 44 | + x = u(2) | ||
| 45 | + end if | ||
| 46 | + first = .not. first | ||
| 47 | + end subroutine | ||
| 48 | + | ||
| 49 | + subroutine randn_vector_n(n, x) | ||
| 50 | + integer, intent(in) :: n | ||
| 51 | + real(dp), intent(out) :: x(n) | ||
| 52 | + integer :: i | ||
| 53 | + do i = 1, size(x) | ||
| 54 | + call randn(x(i)) | ||
| 55 | + end do | ||
| 56 | + end subroutine | ||
| 57 | + | ||
| 58 | + subroutine randn_vector(x) | ||
| 59 | + real(dp), intent(out) :: x(:) | ||
| 60 | + call randn_vector_n(size(x), x) | ||
| 61 | + end subroutine | ||
| 62 | + | ||
| 63 | + subroutine randn_matrix(x) | ||
| 64 | + real(dp), intent(out) :: x(:, :) | ||
| 65 | + call randn_vector_n(size(x), x) | ||
| 66 | + end subroutine | ||
| 67 | + | ||
| 68 | + subroutine rand_gamma0(a, first, fn_val) | ||
| 69 | + ! Returns a psuedorandom scalar drawn from the gamma distribution. | ||
| 70 | + ! | ||
| 71 | + ! The shape parameter a >= 1. | ||
| 72 | + ! | ||
| 73 | + ! [1] Marsaglia, G., & Tsang, W. W. (2000). A Simple Method for Generating | ||
| 74 | + ! Gamma Variables. ACM Transactions on Mathematical Software (TOMS), 26(3), | ||
| 75 | + ! 363–372. | ||
| 76 | + real(dp), intent(in) :: a | ||
| 77 | + logical, intent(in) :: first | ||
| 78 | + real(dp), intent(out) :: fn_val | ||
| 79 | + real(dp), save :: c, d | ||
| 80 | + real(dp) :: U, v, x | ||
| 81 | + if (a < 1) call stop_error("Shape parameter must be >= 1") | ||
| 82 | + if (first) then | ||
| 83 | + d = a - 1._dp/3 | ||
| 84 | + c = 1/sqrt(9*d) | ||
| 85 | + end if | ||
| 86 | + do | ||
| 87 | + do | ||
| 88 | + call randn(x) | ||
| 89 | + v = (1 + c*x)**3 | ||
| 90 | + if (v > 0) exit | ||
| 91 | + end do | ||
| 92 | + call random_number(U) | ||
| 93 | + ! Note: the number 0.0331 below is exact, see [1]. | ||
| 94 | + if (U < 1 - 0.0331_dp*x**4) then | ||
| 95 | + fn_val = d*v | ||
| 96 | + exit | ||
| 97 | + else if (log(U) < x**2/2 + d*(1 - v + log(v))) then | ||
| 98 | + fn_val = d*v | ||
| 99 | + exit | ||
| 100 | + end if | ||
| 101 | + end do | ||
| 102 | + end subroutine | ||
| 103 | + | ||
| 104 | + subroutine rand_gamma_scalar(a, x) | ||
| 105 | + real(dp), intent(in) :: a | ||
| 106 | + real(dp), intent(out) :: x | ||
| 107 | + call rand_gamma0(a, .true., x) | ||
| 108 | + end subroutine | ||
| 109 | + | ||
| 110 | + subroutine rand_gamma_vector_n(a, n, x) | ||
| 111 | + real(dp), intent(in) :: a | ||
| 112 | + integer, intent(in) :: n | ||
| 113 | + real(dp), intent(out) :: x(n) | ||
| 114 | + integer :: i | ||
| 115 | + call rand_gamma0(a, .true., x(1)) | ||
| 116 | + do i = 2, size(x) | ||
| 117 | + call rand_gamma0(a, .false., x(i)) | ||
| 118 | + end do | ||
| 119 | + end subroutine | ||
| 120 | + | ||
| 121 | + subroutine rand_gamma_vector(a, x) | ||
| 122 | + real(dp), intent(in) :: a | ||
| 123 | + real(dp), intent(out) :: x(:) | ||
| 124 | + call rand_gamma_vector_n(a, size(x), x) | ||
| 125 | + end subroutine | ||
| 126 | + | ||
| 127 | + subroutine rand_gamma_matrix(a, x) | ||
| 128 | + real(dp), intent(in) :: a | ||
| 129 | + real(dp), intent(out) :: x(:, :) | ||
| 130 | + call rand_gamma_vector_n(a, size(x), x) | ||
| 131 | + end subroutine | ||
| 132 | + | ||
| 133 | + end module | ||
0 commit comments