Modified random distribution

Given a random number generator, (rng), generating numbers in the range 0.0 .. 1.0 called rgen, for example; and a function modifier(x) taking an number in the same range and generating the probability that the input should be generated, in the same range 0..1; then implement the following algorithm for generating random numbers to the probability given by function modifier:

while True:
    random1 = rgen()
    random2 = rgen()
    if random2 < modifier(random1):
        answer = random1
        break
    endif
endwhile

Task

Create a modifier function that generates a 'V' shaped probability of number generation using something like, for example:

modifier(x) = 2*(0.5 - x) if x < 0.5 else 2*(x - 0.5)

Create a generator of random numbers with probabilities modified by the above function.
Generate >= 10,000 random numbers subject to the probability modification.
Output a textual histogram with from 11 to 21 bins showing the distribution of the random numbers generated.

Show your output here, on this page.

Julia

<lang>using UnicodePlots

modifier(x) = (y = 2x - 1; y < 0 ? -y : y) modrands(rands1, rands2) = [x for (i, x) in enumerate(rands1) if rands2[i] < modifier(x)] histogram(modrands(rand(50000), rand(50000)), nbins = 20)

</lang>

Output:


                ┌                                        ┐ 
   [0.0 , 0.05) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 2412   
   [0.05, 0.1 ) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 2164      
   [0.1 , 0.15) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1850           
   [0.15, 0.2 ) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1652              
   [0.2 , 0.25) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1379                  
   [0.25, 0.3 ) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1121                     
   [0.3 , 0.35) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇ 903                         
   [0.35, 0.4 ) ┤▇▇▇▇▇▇▇▇▇▇ 695                            
   [0.4 , 0.45) ┤▇▇▇▇▇▇ 407                                
   [0.45, 0.5 ) ┤▇▇ 118                                    
   [0.5 , 0.55) ┤▇▇ 126                                    
   [0.55, 0.6 ) ┤▇▇▇▇▇ 358                                 
   [0.6 , 0.65) ┤▇▇▇▇▇▇▇▇▇ 639                             
   [0.65, 0.7 ) ┤▇▇▇▇▇▇▇▇▇▇▇▇ 837                          
   [0.7 , 0.75) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1121                    
   [0.75, 0.8 ) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1332                 
   [0.8 , 0.85) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1608             
   [0.85, 0.9 ) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1920         
   [0.9 , 0.95) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 2204     
   [0.95, 1.0 ) ┤▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 2348   
                └                                        ┘
                                Frequency

Python

<lang python>import random from typing import List, Callable, Optional

def modifier(x: float) -> float:

   """
   V-shaped, modifier(x) goes from 1 at 0 to 0 at 0.5 then back to 1 at 1.0 .

   Parameters
   ----------
   x : float
       Number, 0.0 .. 1.0 .

   Returns
   -------
   float
       Target probability for generating x; between 0 and 1.

   """
   return 2*(.5 - x) if x < 0.5 else 2*(x - .5)

def modified_random_distribution(modifier: Callable[[float], float],

                                n: int) -> List[float]:
   """
   Generate n random numbers between 0 and 1 subject to modifier.

   Parameters
   ----------
   modifier : Callable[[float], float]
       Target random number gen. 0 <= modifier(x) < 1.0 for 0 <= x < 1.0 .
   n : int
       number of random numbers generated.

   Returns
   -------
   List[float]
       n random numbers generated with given probability.

   """
   d: List[float] = []
   while len(d) < n:
       r1 = prob = random.random()
       if random.random() < modifier(prob):
           d.append(r1)
   return d

if __name__ == '__main__':

   from collections import Counter

   data = modified_random_distribution(modifier, 50_000)
   bins = 15
   counts = Counter(d // (1 / bins) for d in data)
   #
   mx = max(counts.values())
   print("   BIN, COUNTS, DELTA: HISTOGRAM\n")
   last: Optional[float] = None
   for b, count in sorted(counts.items()):
       delta = 'N/A' if last is None else str(count - last)
       print(f"  {b / bins:5.2f},  {count:4},  {delta:>4}: "
             f"{'#' * int(40 * count / mx)}")
       last = count</lang>

Output:

   BIN, COUNTS, DELTA: HISTOGRAM

   0.00,  6326,   N/A: ########################################
   0.07,  5327,  -999: #################################
   0.13,  4487,  -840: ############################
   0.20,  3495,  -992: ######################
   0.27,  2601,  -894: ################
   0.33,  1744,  -857: ###########
   0.40,   914,  -830: #####
   0.47,   225,  -689: #
   0.53,   899,   674: #####
   0.60,  1783,   884: ###########
   0.67,  2623,   840: ################
   0.73,  3566,   943: ######################
   0.80,  4383,   817: ###########################
   0.87,  5422,  1039: ##################################
   0.93,  6205,   783: #######################################

REXX

<lang rexx>/*REXX program generates a "V" shaped probability of number generation using a modifier.*/ parse arg randn bins seed . /*obtain optional argument from the CL.*/ if randN== | randN=="," then randN= 100000 /*Not specified? Then use the default.*/ if bins== | bins=="," then bins= 20 /* " " " " " " */ if datatype(seed, 'W') then call random ,,seed /* " " " " " " */ call MRD !.= 0

     do j=1  for randN;   bin= @.j*bins%1
     !.bin= !.bin + 1                           /*bump the applicable bin counter.     */
     end   /*j*/

mx= 0

     do k=1  for randN;   mx= max(mx, !.k)      /*find the maximum, used for histograph*/
     end   /*k*/

say ' bin' say '────── ' center('(with ' commas(randN) " samples", 80 - 10)

      do b=0  for bins;  say format(b/bins,2,2)   copies('■', 70*!.b%mx)" "   commas(!.b)
      end   /*b*/

exit 0 /*──────────────────────────────────────────────────────────────────────────────────────*/ commas: arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ? rand: return random(0, 100000) / 100000 /*──────────────────────────────────────────────────────────────────────────────────────*/ modifier: parse arg y; if y<.5 then return 2 * (.5 - y)

                                 else return  2 * ( y - .5)

/*──────────────────────────────────────────────────────────────────────────────────────*/ MRD: #=0; @.= /*MRD: Modified Random distribution. */

           do until #==randN;      r= rand()    /*generate a random number; assign bkup*/
           if rand()>=modifier(r)  then iterate /*Doesn't meet requirement?  Then skip.*/
           #= # + 1;               @.#= r       /*bump counter; assign the MRD to array*/
           end   /*until*/
         return</lang>

output when using the default inputs:

  bin
──────                         (with  100,000  samples
 0.00 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  9,476
 0.05 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  8,471
 0.10 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  7,528
 0.15 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  6,403
 0.20 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  5,593
 0.25 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  4,541
 0.30 ■■■■■■■■■■■■■■■■■■■■■■■■  3,424
 0.35 ■■■■■■■■■■■■■■■■■■  2,514
 0.40 ■■■■■■■■■■■  1,508
 0.45 ■■■  463
 0.50 ■■■  493
 0.55 ■■■■■■■■■■  1,501
 0.60 ■■■■■■■■■■■■■■■■■■  2,508
 0.65 ■■■■■■■■■■■■■■■■■■■■■■■■  3,416
 0.70 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  4,574
 0.75 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  5,556
 0.80 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  6,506
 0.85 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  7,551
 0.90 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  8,383
 0.95 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  9,590