simulating a trucated normal

[nickg]

I'm interested in simulating a series of random normals truncated at upper and lower limits. In Fortran it is easy to write a code that rejects values beyond the truncation limits using line numbers and IFLG = 0, so that out-of-bounds samples are rejected and a new random is drawn for the same iteration number (and tested again).
Can this be done in C++?
Any code suggestions would be appreciated.

[Repeater]

#include <iostream>
#include <random>
#include <vector>

int main()
{
    // Set up generator
    std::random_device rd{};
    std::mt19937 gen{rd()};
 
    // Normal distribution, mean 10.0, standard dev 3.0
    std::normal_distribution<> d{10.0,3.0};

    // Prepare somewhere to put the numbers
    std::vector<double> randomNumbers;

    // We'll keep going until we have 100 numbers
    while (randomNumbers.size() < 100)
    {
        // generate a random number
	auto randVal = d(gen);

        // Reject numbers below zero and above 20
	if (randVal > 0 && randVal < 20)
	  {
	    randomNumbers.push_back(randVal);
	  }
    }

    // Output the numbers to screen
    for (const auto& val : randomNumbers)
      {
	std::cout << val << ' ';
      }
 }

[lastchance]

This may or may not be what you want (depending on whether you want to realign probabilities to the truncated range).

You have a truly bizarre understanding of how you would do it in modern Fortran.

#include <iostream>
#include <random>
#include <ctime>
using namespace std;

mt19937 gen( time( 0 ) );
normal_distribution<double> Z;     // standard normal Z ~ N(0,1)

//======================================================================

double truncatedNormal( double mean, double sigma, double lower, double upper )
{
   double X = mean + sigma * Z( gen );      // Rescale standard normal to desired parameters
   if ( X >= lower && X <= upper )          // If in range, return it
   {
      return X;
   }
   else                                     // Otherwise, recursively try again
   {
      return truncatedNormal( mean, sigma, lower, upper );
   }
}

//======================================================================

int main()
{
   const int N = 10;                       // Number of samples
   const double mean = 100, sigma = 10;    // mean and standard deviation of un-truncated distribution
   const double lower = 95, upper = 115;   // Cut-offs

   for ( int i = 0; i < N; i++ ) cout << truncatedNormal( mean, sigma, lower, upper ) << '\n';
}

[nickg]

Thanks. Please illustrate how you would cod eit in modern Fortran

[lastchance]

nickg wrote:

simulating a series of random normals truncated at upper and lower limits. In Fortran it is easy to write a code that rejects values beyond the truncation limits using line numbers and IFLG = 0

nickg wrote:

Thanks. Please illustrate how you would code it in modern Fortran

module Distribution
   implicit none
   real, parameter :: PI = 3.14159

contains
   subroutine standardNormal( Z )                ! Uses Box-Muller method
      real, intent(out) :: Z
      real u(2)

      call random_number( u )                    ! 2 variables uniformly distributed on (0,1)
      u(1) = max( u(1), 1.0e-30 )                ! paranoia
      Z = sqrt( -2.0 * log( u(1) ) ) * cos( 2.0 * PI * u(2) )
   end subroutine standardNormal

   real function truncatedNormal( mean, sigma, lower, upper )
      real, intent(in) :: mean, sigma, lower, upper
      real Z

      do
         call standardNormal( Z )
         truncatedNormal = mean + sigma * Z
         if ( truncatedNormal >= lower .and. truncatedNormal <= upper ) return
      end do
   end function truncatedNormal

end module Distribution


program main
   use Distribution                              ! access module
   implicit none

   integer :: N = 10
   real :: mean = 100, sigma = 10, lower = 95, upper = 115
   integer i

   call random_seed
   write( *, "( f6.2 )" ) ( truncatedNormal( mean, sigma, lower, upper ), i = 1, N )

end program main

https://onlinegdb.com/Bk9To7zBU

[nickg]

Thanks you once again