Loading ...
Sorry, an error occurred while loading the content.
 

Generate random numbers following a normal (Gaussian) distribution

Expand Messages
  • citromatik
    Hi all, I am trying to get random numbers that follows a Gaussian distribution given the mean and standard deviation of the distribution. So far, I haven t
    Message 1 of 6 , Sep 3, 2009
      Hi all,

      I am trying to get random numbers that follows a Gaussian distribution given
      the mean and standard deviation of the distribution.

      So far, I haven't found any OCaml code to do this, so I translated a Perl
      module that does exactly what I wanted:
      http://search.cpan.org/~dagolden/Math-Random-OO-0.21/lib/Math/Random/OO/Normal.pm

      The algorithm used in the module can be found here:
      http://home.online.no/~pjacklam/notes/invnorm/

      This is the OCaml code I wrote, in case anyone is interested:

      (* This code generates random normals from the inverse of the cumulative
      normal distribution using an approximation algorithm developed by Peter J.
      Acklam and released into the public domain. This algorithm claims a relative
      error of less than 1.15e-9 over the entire region. *)
      (* See http://home.online.no/~pjacklam/notes/invnorm/ for details *)

      exception OutOfRange

      (* Coefficients in rational approximations. *)
      let a = [|-3.969683028665376e+01;2.209460984245205e+02;
      -2.759285104469687e+02;1.383577518672690e+02;
      -3.066479806614716e+01;2.506628277459239e+00|]

      let b = [|-5.447609879822406e+01;1.615858368580409e+02;
      -1.556989798598866e+02;6.680131188771972e+01;
      -1.328068155288572e+01|]

      let c = [|-7.784894002430293e-03;-3.223964580411365e-01;
      -2.400758277161838e+00;-2.549732539343734e+00;
      4.374664141464968e+00;2.938163982698783e+00|]

      let d = [|7.784695709041462e-03;3.224671290700398e-01;
      2.445134137142996e+00;3.754408661907416e+00|]

      (* Define break-points *)
      let p_low = 0.02425
      let p_high = 1.0 -. p_low

      let lqtnorm = function
      | p when p > 0.0 && p < p_low ->
      let q = sqrt (-2.0 *. log p) in
      (((((c.(0)*.q+.c.(1))*.q+.c.(2))*.q+.c.(3))*.q+.c.(4))*.q+.c.(5)) /.
      ((((d.(0)*.q+.d.(1))*.q+.d.(2))*.q+.d.(3))*.q+.1.0)
      | p when p >= p_low && p <= p_high ->
      let q = p -. 0.5 in
      let r = q *. q in
      (((((a.(0)*.r+.a.(1))*.r+.a.(2))*.r+.a.(3))*.r+.a.(4))*.r+.a.(5))*.q
      /.
      (((((b.(0)*.r+.b.(1))*.r+.b.(2))*.r+.b.(3))*.r+.b.(4))*.r+.1.0)
      | p when p > p_high && p < 1.0 ->
      let q = sqrt (-2.0 *. log (1.0 -. p)) in
      -.(((((c.(0)*.q+.c.(1))*.q+.c.(2))*.q+.c.(3))*.q+.c.(4))*.q+.c.(5)) /.
      ((((d.(0)*.q+.d.(1))*.q+.d.(2))*.q+.d.(3))*.q+.1.0)
      | _ -> raise OutOfRange

      class gaussian (mean:float) (stdev:float) =
      object
      initializer Random.self_init ()
      method set_seed i = Random.init i
      method next =
      let stdev = abs_float stdev in
      let rnd = let n' = Random.float 1.0 in if n' = 0.0 then 1e-254 else
      n' in
      let norm =
      if rnd <= 0.5 then lqtnorm(rnd)
      else -.(lqtnorm(1.0 -. rnd)) in
      norm *. stdev +. mean
      end

      Does anybody know of any other OCaml module for doing this kind of things?
      Do you think that this code can be submitted to any OCaml code repository?

      Thanks in advance,

      M;

      --
      View this message in context: http://www.nabble.com/Generate-random-numbers-following-a-normal-%28Gaussian%29-distribution-tp25275888p25275888.html
      Sent from the Ocaml Beginner mailing list archive at Nabble.com.
    • Hezekiah M. Carty
      ... ... The GNU Science Library (GSL) has several random number distribution libraries, including Gaussian. There are OCaml bindings for GSL available
      Message 2 of 6 , Sep 3, 2009
        On Thu, Sep 3, 2009 at 8:20 AM, citromatik<miguel.pignatelli@...> wrote:
        > I am trying to get random numbers that follows a Gaussian distribution given
        > the mean and standard deviation of the distribution.
        >
        > So far, I haven't found any OCaml code to do this, so I translated a Perl
        > module that does exactly what I wanted:
        > http://search.cpan.org/~dagolden/Math-Random-OO-0.21/lib/Math/Random/OO/Normal.pm
        >
        > The algorithm used in the module can be found here:
        > http://home.online.no/~pjacklam/notes/invnorm/
        >
        <snip>
        >
        > Does anybody know of any other OCaml module for doing this kind of things?
        > Do you think that this code can be submitted to any OCaml code repository?

        The GNU Science Library (GSL) has several random number distribution
        libraries, including Gaussian. There are OCaml bindings for GSL
        available here:

        http://code.google.com/p/ocamlgsl/

        Packages for the GSL bindings are included in Debian, Ubuntu and
        Fedora (and possibly others).

        You can submit a code snippet on forge.ocamlcore.org:

        http://forge.ocamlcore.org/snippet/browse.php?by=lang&lang=1

        Hez
      • Richard Jones
        ... Note that perl4caml lets you use Perl modules directly inside OCaml code. Of course it won t be very fast, since it still interprets the Perl, but might
        Message 3 of 6 , Sep 3, 2009
          On Thu, Sep 03, 2009 at 06:20:53AM -0700, citromatik wrote:
          > So far, I haven't found any OCaml code to do this, so I translated a Perl
          > module that does exactly what I wanted:
          > http://search.cpan.org/~dagolden/Math-Random-OO-0.21/lib/Math/Random/OO/Normal.pm

          Note that perl4caml lets you use Perl modules directly inside OCaml
          code. Of course it won't be very fast, since it still interprets the
          Perl, but might be more convenient than rewriting the code.

          http://merjis.com/developers/perl4caml

          perl4caml is packaged in the major Linux distros.

          Rich.

          --
          Richard Jones
          Red Hat
        • Richard Jones
          ... OCamlForge is a hosting service you can use for new projects: http://forge.ocamlcore.org/ Alternately, just host your code on the usual places (gitorious,
          Message 4 of 6 , Sep 3, 2009
            On Thu, Sep 03, 2009 at 06:20:53AM -0700, citromatik wrote:
            > Do you think that this code can be submitted to any OCaml code repository?

            OCamlForge is a hosting service you can use for new projects:

            http://forge.ocamlcore.org/

            Alternately, just host your code on the usual places (gitorious,
            Google Code, berlios, etc):

            http://en.wikipedia.org/wiki/Category:Free_software_hosting_facilities

            Rich.

            --
            Richard Jones
            Red Hat
          • citromatik
            Hi Richard, Thanks for the tip, I know perl4caml (very cool indeed!), but in this case efficiency is needed and the translation to OCaml is fairly easy. Thanks
            Message 5 of 6 , Sep 3, 2009
              Hi Richard,

              Thanks for the tip, I know perl4caml (very cool indeed!), but in this case
              efficiency is needed and the translation to OCaml is fairly easy.

              Thanks anyway!

              M;




              Richard Jones-4 wrote:
              >
              > On Thu, Sep 03, 2009 at 06:20:53AM -0700, citromatik wrote:
              >> So far, I haven't found any OCaml code to do this, so I translated a Perl
              >> module that does exactly what I wanted:
              >> http://search.cpan.org/~dagolden/Math-Random-OO-0.21/lib/Math/Random/OO/Normal.pm
              >
              > Note that perl4caml lets you use Perl modules directly inside OCaml
              > code. Of course it won't be very fast, since it still interprets the
              > Perl, but might be more convenient than rewriting the code.
              >
              > http://merjis.com/developers/perl4caml
              >
              > perl4caml is packaged in the major Linux distros.
              >
              > Rich.
              >
              > --
              > Richard Jones
              > Red Hat
              >
              >

              --
              View this message in context: http://www.nabble.com/Generate-random-numbers-following-a-normal-%28Gaussian%29-distribution-tp25275888p25284573.html
              Sent from the Ocaml Beginner mailing list archive at Nabble.com.
            • Yang Ye
              Box-muller method is a common way to get random number of normal distribution. A sample of code is here (taken from
              Message 6 of 6 , Sep 4, 2009
                Box-muller method is a common way to get random number of normal
                distribution. A sample of code is here (taken from
                http://www.uncarved.com/blog/practical_ocaml.mrk)

                open Random;;

                (* initialize the random number generator *)
                Random.self_init;;

                (* get a random gaussian using a Box-Muller transform, described
                * here http://en.wikipedia.org/wiki/Box-Muller_transform *)
                let rec get_one_gaussian_by_box_muller () =
                (* Generate two uniform numbers from -1 to 1 *)
                let x = Random.float 2.0 -. 1.0 in
                let y = Random.float 2.0 -. 1.0 in
                let s = x*.x +. y*.y in
                if s > 1.0 then get_one_gaussian_by_box_muller ()
                else x *. sqrt (-2.0 *. (log s) /. s)
                ;;


                Regards,
                Yang Ye

                On Fri, Sep 4, 2009 at 5:40 AM, citromatik <miguel.pignatelli@...> wrote:

                >
                >
                >
                > Hi Richard,
                >
                > Thanks for the tip, I know perl4caml (very cool indeed!), but in this case
                > efficiency is needed and the translation to OCaml is fairly easy.
                >
                > Thanks anyway!
                >
                > M;
                >
                >
                > Richard Jones-4 wrote:
                > >
                > > On Thu, Sep 03, 2009 at 06:20:53AM -0700, citromatik wrote:
                > >> So far, I haven't found any OCaml code to do this, so I translated a
                > Perl
                > >> module that does exactly what I wanted:
                > >>
                > http://search.cpan.org/~dagolden/Math-Random-OO-0.21/lib/Math/Random/OO/Normal.pm<http://search.cpan.org/%7Edagolden/Math-Random-OO-0.21/lib/Math/Random/OO/Normal.pm>
                > >
                > > Note that perl4caml lets you use Perl modules directly inside OCaml
                > > code. Of course it won't be very fast, since it still interprets the
                > > Perl, but might be more convenient than rewriting the code.
                > >
                > > http://merjis.com/developers/perl4caml
                > >
                > > perl4caml is packaged in the major Linux distros.
                > >
                > > Rich.
                > >
                > > --
                > > Richard Jones
                > > Red Hat
                > >
                > >
                >
                > --
                > View this message in context:
                > http://www.nabble.com/Generate-random-numbers-following-a-normal-%28Gaussian%29-distribution-tp25275888p25284573.html
                > Sent from the Ocaml Beginner mailing list archive at Nabble.com.
                >
                >
                >



                --
                Regards,
                Yang Ye


                [Non-text portions of this message have been removed]
              Your message has been successfully submitted and would be delivered to recipients shortly.