CODEWITHSUNDEEP
fortranrandomnormal-distribution
ryanjdillon
ryanjdillon

Reputation: 18978

Generate normally distributed data within a range using Fortran

I would like to generate an array of normally distributed data around a given mean and withing given upper and lower bounds.

I found a function here that uses the Marsaglia polar method (a.k.a polar form of the Box-Müller transformation).

How might I modify this to generate points within -10 and +10 integer values around the mean?

  function normal(mean, sigma)

  ! mean  : mean of distribution
  ! sigma : number of standard deviations

  implicit none

  integer, parameter:: b8 = selected_real_kind(14)
  real(b8), parameter :: pi = 3.141592653589793239_b8
  real(b8) normal
  real rand_num
  real(b8) tmp
  real(b8) mean
  real(b8) sigma
  real(b8) fac
  real(b8) gsave
  real(b8) rsq
  real(b8) r1
  real(b8) r2
  integer flag
  save flag
  save gsave
  data flag /0/

  if (flag.eq.0) then
    rsq=2.0_b8

    do while(rsq.ge.1.0_b8.or.rsq.eq.0.0_b8) ! new from for do
      call random_number(rand_num)
      r1=2.0_b8*rand_num-1.0_b8
      call random_number(rand_num)
      r2=2.0_b8*rand_num-1.0_b8
      rsq=r1*r1+r2*r2
    enddo

    fac=sqrt(-2.0_b8*log(rsq)/rsq)
    gsave=r1*fac
    tmp=r2*fac
    flag=1
  else
    tmp=gsave
    flag=0
  endif

  normal=tmp*sigma+mean

  return
  endfunction normal

Upvotes: 1

Views: 3207

Answers (1)

francescalus
francescalus

Reputation: 32451

There are different ways to sample from a truncated normal distribution. The "best" way depends on your mean and variance.

A simple way is just plain rejection: generate a deviate and if it's too big or small throw it away and repeat. If your parameters are such that you don't reject many this is a good method.

For the routine specified, doing this rejection is simple: sticking the whole calculation part (including the flag test) into a do...end do loop with an if (abs(normal).lt.10) exit should do it.

[This isn't elegant, and extra bookkeeping may be required, with samples being generated as pairs. But, again, if rejection happens rarely that's not too bad. Tidying is left as an exercise for the reader.]

That's how one may change the routine, but as @george comments, that probably isn't the best approach, even with the rejection. Certainly, having a function called normal which takes mean and sigma and returns a deviate which isn't from a normal distribution could be confusing.

function truncated_normal(mu, sigma, lim_a, lim_b) !Very side-effecty
  ... declarations ...
  if ( ... test for wanting the rejection method ... ) then
    do
      truncated_normal = normal(mu, sigma)
      if (truncated_normal.gt.lim_a.and.truncated_normal.lt.lim_b) exit
    end do
  else
    ... the other magical method when rejection is too likely ...
  end if
end

Upvotes: 3

Related Questions