4
$\begingroup$

I am working on a Fortran code that uses a uniformly spaced grid in two directions (x,y). Which works fine, but when I need to study a certain problem with good resolution, I need to increase the number of points in one of the directions. This works fine, but it requires a lot of computational resources.

I want to implement a non-uniform grid in one of those directions (y direction), in particular, a Gaussian spacing around a mean value. In other words, I need a bunch of points around a certain y value, then points with an increasingly larger space as we go further from this point (for both lower and higher values of y).

I have looked around but can't seem to find anything useful. Maybe I am using the wrong keywords. I found many things for Python or Matlab that use the already implemented logspace or other functions. I also found a logspace in Fortran (I did not find the source/algorithm), but I would prefer to code my own. If it has already been solved, please let me know.

$\endgroup$

1 Answer 1

9
$\begingroup$

What you need is the Gaussian quantile function. The quantile function outputs the value of a random variable such that its probability is less than or equal to an input probability value.

The quantile function of the standard normal distribution is called the probit function, and can be expressed in terms of the inverse error function:

$$ \Phi^{-1}(p) = {\sqrt{2}}\texttt{erf}^{-1}(2p - 1), \quad p \in (0,1) ~. $$

For a normal random variable with mean $\mu$ and variance $\sigma^{2}$, the quantile function is,

$$ F^{-1}(p) = \mu + \sigma \Phi^{-1}(p) = \mu + \sigma \sqrt{2} \texttt{erf}^{-1}(2p - 1), \quad p \in (0,1) ~. $$

So, for example, to generate a set of 100 Gaussian-spaced points with mean $\mu$ and standard deviation $\sigma$, you can pass a set of 100 uniformly-spaced values $p$ in the range $(0,1)$ to the Gaussian quantile function above. Implementing the Gaussian quantile function is a rather involved task because it requires implementing the inverse error function.

The ParaMonte Fortran library has generic (arbitrary-precision) implementations of both Gaussian quantile function and the inverse error function. Both functionalities are documented in getNormQuan and getErfInv.

If you need linear/log-spaced vector, check out the pm_arraySpace module of the ParaMonte library. Again, example usage and documentation can be found in pm_arraySpace.

p.s. I am the original developer of the ParaMonte library. If you need help extracting these functionalities for use in your codebase, please open an issue in the project's GitHub repository to further assist you.

$\endgroup$
4
  • $\begingroup$ Excellent Ill check that. The code I use already has a function for the error function (which is not hard to implement). I'll go check the documentation you gave me, and use it if I find it too hard to implement one by myself. Thank you. $\endgroup$
    – Gundro
    Commented Jun 19 at 13:56
  • $\begingroup$ To be clear, you need the inverse of the error function. Incidentally, the Fortran language has the intrinsic erf() Error Function, so you don't need an external implementation but it does not have the inverse. $\endgroup$
    – Scientist
    Commented Jun 19 at 14:01
  • $\begingroup$ Yes I understood, it is just that the original developer implemented a different version of the erf which from what I understood is better in the regime we work in, compared tot he intrinsic erf. I am checking the inverse rightnow. $\endgroup$
    – Gundro
    Commented Jun 19 at 14:57
  • $\begingroup$ I want to thank you. I implemented an inv_erf on my own and then used it as you mentioned. It works so well that I now have a Gaussianly distributed array instead of a uniform one. Thank you very much. I just need to implement it in the main code, and it should work after a few tests. $\endgroup$
    – Gundro
    Commented Jun 20 at 13:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.