Sampling from an arbitrary continuous distribution

1

Is there a Fortran package that will generate random variates from an arbitrary probability distribution defined by a function, for example

function f(x) result(y)
real, intent(in) :: x
real             :: y
y = 1/(2*exp(abs(x))) ! Laplace density
end function f

There are codes to generate from common distributions such as the normal, Student’s t, and Laplace, but I would like to be able to simulate from any function that is always nonnegative and integrates to 1 from -Inf to Inf. There is the method of rejection sampling, but how do you automatically choose a proposal distribution? Another approach is to compute the cumulative distribution function (CDF) at many points by integrating function f and then obtaining the inverse CDF by root-finding, as discussed here.

1 Like
2

That is exactly the reason why I (and other team members) developed the ParaMonte Fortran library; for sampling positive functions of arbitrary shape in arbitrary dimensions. Here are the steps:

  1. Download the ParaMonte prebuilt library for Windows, Linux, or macOS, or build it from scratch if needed for Windows, Linux, or macOS.

  2. Relpace the contents of the logfunc.f90 file with your objective function implementation, like the Laplace density:

module LogFunc_mod
    use paramonte, only: IK, RK
    implicit none
    integer(IK), parameter :: NDIM = 1_IK
contains
    !> \brief
    !> Return the natural logarithm of the Laplace density: 1._RK / (2_IK * exp(abs(Point(1)))).
    function getLogFunc(ndim,Point) result(logFunc)
        integer(IK), intent(in) :: ndim
        real(RK), intent(in)    :: Point(ndim)
        real(RK), parameter     :: NEG_LOG_TWO = -log(2._RK)
        real(RK)                :: logFunc
        logFunc = NEG_LOG_TWO - abs(Point(1))
    end function getLogFunc
end module LogFunc_mod

  1. Call the build.bat (on Windows) or build.sh && ./run.sh (on Linux / macOS) to build and run the sampler. This will generate a set of files one of which (*_sample.txt) contains the final set of random samples from the objective function.

  2. We do not yet have a flexible visualization library in Fortran. But the ParaMonte library has dedicated Python and MATLAB interfaces that can readily post-process the results and create professional research-quality visualizations from the output files with minimal coding. For example in Jupyter Python,

!pip install --user paramonte
import paramonte as pm
pmpd = pm.Paradram() # generate an instance of the ParaMonte::ParaDRAM sampler.
sample = pmpd.readSample("./out/")[0] # read the output sample. The program is smart to find the file in this folder.
sample.plot.histplot() # generate a histogram plot.
sample.plot.line() # generate a line plot.

The above code will generate two plots like the following,

laplace_hist
laplace_hist876×576 11.6 KB

laplace_trace
laplace_trace870×576 73.7 KB

There is a lot more that can be done with this library, both in terms of sampling and in terms of visualization.

I originally created this library for my own research needs, but decided to polish and release it for public use by the Fortran community (and Python, MATLAB, C, C++, R, …). I’d be very happy to help you further with implementing your objective functions and you end up using it, your feedbacks are also appreciated very much.

7 Likes
3

An absolutely interesting project @shahmoradi, congratulations! Is there any plans to make ParaMonte available through fpm? Are there any fpm missing features that are necessary to make this happen?

4

Thank you @epagone. Automatic build with fpm is definitely on the priority todo list. I have not tested fpm yet so I am not sure what incompatibilities may exist. One thing I know for sure is that the module naming convention of the library should be already mostly consistent with fpm. Aside from fpm, we have also developed fully automated build scripts for both Windows and Unix systems that do not leave any manual work for the user, even the installation of basic requirements like Fortran compiler, CMake, MPI, and (Open)Coarray(s) Fortran libraries are automatically checked and done as part of the ParaMonte library build. Such a level of automation is really needed for higher-level users of the library who typically do not even know what Fortran is, let alone the installation of prereqs.
And one last thing, although I do not have the official benchmarks yet, internal testings make me confident this is library is among the fastest (if not the fastest), most comprehensive, and most reliable adaptive MCMC samplers ever developed in any programming language. So use it with confidence and peace of mind :-).

2 Likes