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