# Sampling from an arbitrary continuous distribution

**URL:** <https://fortran-lang.discourse.group/t/sampling-from-an-arbitrary-continuous-distribution/1381>\
**Category:** Uncategorized\
**Created:** [June 11, 2021, 11:44pm UTC](https://fortran-lang.discourse.group/t/sampling-from-an-arbitrary-continuous-distribution/1381 "2021-06-11T23:44:59Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [June 11, 2021, 11:44pm UTC](https://fortran-lang.discourse.group/t/sampling-from-an-arbitrary-continuous-distribution/1381/1 "2021-06-11T23:44:59Z")

</div>

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

```auto
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](https://en.wikipedia.org/wiki/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](https://blogs.sas.com/content/iml/2013/07/22/the-inverse-cdf-method.html).

---

<div class="post-metadata">

**Author:** ![shahmoradi](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/shahmoradi/32/3151_2.png) [@shahmoradi](https://fortran-lang.discourse.group/u/shahmoradi)\
**Post date:** [June 12, 2021, 2:41am UTC](https://fortran-lang.discourse.group/t/sampling-from-an-arbitrary-continuous-distribution/1381/2 "2021-06-12T02:41:33Z")

</div>

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

1. Download the ParaMonte prebuilt library for [Windows](https://www.cdslab.org/paramonte/notes/installation/windows/#using-the-paramonte-library-prebuilt-via-the-intel-compiler-suite), [Linux](https://www.cdslab.org/paramonte/notes/installation/linux/#using-the-paramonte-library-prebuilt-via-the-intel-compiler-suite), or [macOS](https://www.cdslab.org/paramonte/notes/installation/macos/#using-the-paramonte-library-prebuilt-via-the-intel-compiler-suite), or build it from scratch if needed for [Windows](https://www.cdslab.org/paramonte/notes/installation/windows/#building-the-paramonte-library-locally-on-your-system), [Linux](https://www.cdslab.org/paramonte/notes/installation/linux/#building-the-paramonte-library-locally-on-your-system), or [macOS](https://www.cdslab.org/paramonte/notes/installation/windows/#building-the-paramonte-library-locally-on-your-system).

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

```auto
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](https://www.cdslab.org/paramonte/notes/installation/python/) and [MATLAB](https://www.cdslab.org/paramonte/notes/installation/matlab/) interfaces that can readily post-process the results and create [professional research-quality visualizations](https://www.cdslab.org/paramonte/notes/examples/matlab/mlx/sampling_multivariate_normal_distribution_via_paradram.html) from the output files with minimal coding. For example in Jupyter Python,

```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](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/1X/d78e2945797be8dd6ba503895946e5acf0b51616.png)  
 ![laplace_trace](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/1X/e5a6fd207392675bbf4def82858abc099bee5009.png)

There is a lot more that can be done with this library, both [in terms of sampling](https://www.cdslab.org/paramonte/notes/examples/python/jupyter/) and [in terms of visualization](https://www.cdslab.org/paramonte/notes/examples/matlab/mlx/).

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.

---

<div class="post-metadata">

**Author:** ![epagone](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/epagone/32/252_2.png) [@epagone](https://fortran-lang.discourse.group/u/epagone)\
**Post date:** [June 12, 2021, 10:29am UTC](https://fortran-lang.discourse.group/t/sampling-from-an-arbitrary-continuous-distribution/1381/3 "2021-06-12T10:29:45Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![shahmoradi](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/shahmoradi/32/3151_2.png) [@shahmoradi](https://fortran-lang.discourse.group/u/shahmoradi)\
**Post date:** [June 12, 2021, 11:10am UTC](https://fortran-lang.discourse.group/t/sampling-from-an-arbitrary-continuous-distribution/1381/4 "2021-06-12T11:10:50Z")

</div>

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 :-).