Add rejection sampling

examples/CMakeLists.txt

@@ -42,4 +42,10 @@ target_link_libraries(bootstrap_regression_example fstats)
 add_executable(box_muller_example box_muller_example.f90)
 target_link_libraries(box_muller_example fstats)
 target_link_libraries(box_muller_example ${fplot_LIBRARY})
-target_include_directories(box_muller_example PUBLIC ${fplot_INCLUDE_DIR})
+target_include_directories(box_muller_example PUBLIC ${fplot_INCLUDE_DIR})
+
+# Rejection Sampling Example
+add_executable(rejection_sample_example rejection_sample_example.f90)
+target_link_libraries(rejection_sample_example fstats)
+target_link_libraries(rejection_sample_example ${fplot_LIBRARY})
+target_include_directories(rejection_sample_example PUBLIC ${fplot_INCLUDE_DIR})

examples/rejection_sample_example.f90

@@ -0,0 +1,34 @@
+program example
+    use iso_fortran_env
+    use fstats
+    use fplot_core
+    implicit none
+
+    ! Description:
+    ! This example uses rejection sampling to sample an F distribution.
+
+    ! Local Variables
+    integer(int32), parameter :: npts = 10000
+    real(real64), parameter :: d1 = 1.0d2
+    real(real64), parameter :: d2 = 1.0d2
+    real(real64), parameter :: xmin = 0.0d0
+    real(real64), parameter :: xmax = 2.0d0
+    real(real64) :: x(npts)
+    type(f_distribution) :: dist
+
+    ! Plot Variables
+    type(plot_2d) :: plt
+    type(plot_data_histogram) :: pd
+
+    ! Perform the sampling
+    dist%d1 = d1
+    dist%d2 = d2
+    x = rejection_sample(dist, npts, xmin, xmax)
+
+    ! Plot the resulting distribution
+    call plt%initialize()
+    call pd%define_data(x)
+    call pd%set_transparency(0.5)
+    call plt%push(pd)
+    call plt%draw()
+end program

src/fstats.f90

@@ -55,6 +55,7 @@ module fstats
     public :: bootstrap_regression_statistics
     public :: bootstrap_linear_least_squares
     public :: box_muller_sample
+    public :: rejection_sample
     public :: FS_LEVENBERG_MARQUARDT_UPDATE
     public :: FS_QUADRATIC_UPDATE
     public :: FS_NIELSEN_UPDATE
@@ -2204,17 +2205,6 @@ module function box_muller_sample_real64(mu, sigma) result(rst)
                 !! The pair of random values.
         end function
 
-        module function box_muller_sample_real32(mu, sigma) result(rst)
-            !! Generates a pair of independent, standard, normally distributed
-            !! random values using the Box-Muller transform.
-            real(real32), intent(in) :: mu
-                !! The mean of the distribution.
-            real(real32), intent(in) :: sigma
-                !! The standard deviation of the distribution.
-            real(real32) :: rst(2)
-                !! The pair of random values.
-        end function
-
         module function box_muller_array_real64(mu, sigma, n) result(rst)
             !! Generates an array of normally distributed random values sampled
             !! by the Box-Muller transform.
@@ -2228,26 +2218,26 @@ module function box_muller_array_real64(mu, sigma, n) result(rst)
                 !! A 2N-element array containing the N Box-Muller pairs.
         end function
 
-        module function box_muller_array_real32(mu, sigma, n) result(rst)
-            !! Generates an array of normally distributed random values sampled
-            !! by the Box-Muller transform.
-            real(real32), intent(in) :: mu
-                !! The mean of the distribution.
-            real(real32), intent(in) :: sigma
-                !! The standard deviation of the distribution.
+        module function rejection_sample(tdist, n, xmin, xmax) result(rst)
+            !! Uses rejection sampling to randomly sample a target distribution.
+            class(distribution), intent(in) :: tdist
+                !! The distribution to sample
             integer(int32), intent(in) :: n
-                !! The number of Box-Muller pairs to generate.
-            real(real32), allocatable, dimension(:) :: rst
-                !! A 2N-element array containing the N Box-Muller pairs.
+                !! The number of samples to make.
+            real(real64), intent(in) :: xmin
+                !! The minimum range to explore.
+            real(real64), intent(in) :: xmax
+                !! The maximum range to explore.
+            real(real64), allocatable, dimension(:) :: rst
+                !! An N-element array containing the N samples from the 
+                !! distribution.
         end function
     end interface
 
     interface box_muller_sample
         !! Generates random, normally distributed values via the Box-Muller 
         !! transform.
         module procedure :: box_muller_sample_real64
-        module procedure :: box_muller_sample_real32
         module procedure :: box_muller_array_real64
-        module procedure :: box_muller_array_real32
     end interface
 end module

src/sampling.f90

@@ -28,27 +28,6 @@ module function box_muller_sample_real64(mu, sigma) result(rst)
     rst = [real(z, real64), aimag(z)]
 end function
 
-! --------------------
-module function box_muller_sample_real32(mu, sigma) result(rst)
-    ! Arguments
-    real(real32), intent(in) :: mu
-    real(real32), intent(in) :: sigma
-    real(real32) :: rst(2)
-
-    ! Parameters
-    complex(real32), parameter :: j = (0.0, 1.0)
-
-    ! Local Variables
-    real(real32) :: u1, u2
-    complex(real32) :: z
-
-    ! Process
-    call random_number(u1)
-    call random_number(u2)
-    z = sqrt(-log(u1)) * exp(j * twopi_f * u2)
-    rst = [real(z, real32), aimag(z)]
-end function
-
 ! ------------------------------------------------------------------------------
 module function box_muller_array_real64(mu, sigma, n) result(rst)
     ! Arguments
@@ -71,31 +50,59 @@ module function box_muller_array_real64(mu, sigma, n) result(rst)
     end do
 end function
 
-! --------------------
-module function box_muller_array_real32(mu, sigma, n) result(rst)
+! ******************************************************************************
+! REJECTION SAMPLING
+! ------------------------------------------------------------------------------
+module function rejection_sample(tdist, n, xmin, xmax) result(rst)
     ! Arguments
-    real(real32), intent(in) :: mu
-    real(real32), intent(in) :: sigma
+    class(distribution), intent(in) :: tdist
     integer(int32), intent(in) :: n
-    real(real32), allocatable, dimension(:) :: rst
+    real(real64), intent(in) :: xmin, xmax
+    real(real64), allocatable, dimension(:) :: rst
+
+    ! Parameters
+    real(real64), parameter :: zero = 0.0d0
+    real(real64), parameter :: c_start = 1.01d0
 
     ! Local Variables
-    integer(int32) :: i
+    integer(int32) :: i, j, jmax
+    real(real64) :: u, c, g, f, rng
 
     ! Process
-    if (n < 1) then
-        allocate(rst(0))
-        return
-    end if
-    allocate(rst(2 * n))
-    do i = 1, n
-        rst(2*i-1:2*i) = box_muller_sample(mu, sigma)
+    i = 0
+    j = 0
+    jmax = 1000 * n  ! Guard against insanity
+    rng = xmax - xmin
+    c = c_start
+    allocate(rst(n), source = zero)
+    do while (i <= n)
+        ! Update the acceptance threshold
+        call random_number(u)
+
+        ! Sample from the proposal distribution
+        call random_number(g)
+        g = g * rng + xmin
+
+        ! Sample the target distribution
+        f = tdist%pdf(g)
+
+        ! Test
+        if (u <= f / (c * g)) then
+            i = i + 1
+            rst(i) = g
+        end if
+
+        ! Update C
+        c = max(c, f / g)
+
+        ! Update the infinite loop guard variable
+        j = j + 1
+        if (j == jmax) exit
     end do
 end function
 
-
 ! ******************************************************************************
-! MAXIMUM ENTROPY
+! MAXIMUM ENTROPY SAMPLING
 ! ------------------------------------------------------------------------------