Modified random distribution
You are encouraged to solve this task according to the task description, using any language you may know.
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.
Ada
<lang Ada>with Ada.Text_Io; with Ada.Numerics.Float_Random; with Ada.Strings.Fixed;
procedure Modified_Distribution is
Observations : constant := 20_000; Buckets : constant := 25; Divider : constant := 12; Char : constant Character := '*';
generic with function Modifier (X : Float) return Float; package Generic_Random is function Random return Float; end Generic_Random;
package body Generic_Random is package Float_Random renames Ada.Numerics.Float_Random; Generator : Float_Random.Generator;
function Random return Float is Random_1 : Float; Random_2 : Float; begin loop Random_1 := Float_Random.Random (Generator); Random_2 := Float_Random.Random (Generator); if Random_2 < Modifier (Random_1) then return Random_1; end if; end loop; end Random;
begin Float_Random.Reset (Generator); end Generic_Random;
generic Buckets : in Positive; package Histograms is type Bucket_Index is new Positive range 1 .. Buckets; Bucket_Width : constant Float := 1.0 / Float (Buckets); procedure Clean; procedure Increment_Bucket (Observation : Float); function Observations_In (Bucket : Bucket_Index) return Natural; function To_Bucket (X : Float) return Bucket_Index; function Range_Image (Bucket : Bucket_Index) return String; end Histograms;
package body Histograms is Hist : array (Bucket_Index) of Natural := (others => 0);
procedure Clean is begin Hist := (others => 0); end Clean;
procedure Increment_Bucket (Observation : Float) is Bin : constant Bucket_Index := To_Bucket (Observation); begin Hist (Bin) := Hist (Bin) + 1; end Increment_Bucket;
function Observations_In (Bucket : Bucket_Index) return Natural is (Hist (Bucket));
function To_Bucket (X : Float) return Bucket_Index is (1 + Bucket_Index'Base (Float'Floor (X / Bucket_Width)));
function Range_Image (Bucket : Bucket_Index) return String is package Float_Io is new Ada.Text_Io.Float_Io (Float); Image : String := "F.FF..L.LL"; First : constant Float := Float (Bucket - 1) / Float (Buckets); Last : constant Float := Float (Bucket - 1 + 1) / Float (Buckets); begin Float_Io.Put (Image (1 .. 4), First, Exp => 0, Aft => 2); Float_Io.Put (Image (7 .. 10), Last, Exp => 0, Aft => 2); return Image; end Range_Image;
begin Clean; end Histograms;
function Modifier (X : Float) return Float is (if X in Float'First .. 0.5 then 2.0 * (0.5 - X) else 2.0 * (X - 0.5));
package Modified_Random is new Generic_Random (Modifier => Modifier);
package Histogram_20 is new Histograms (Buckets => Buckets);
function Column (Height : Natural; Char : Character) return String renames Ada.Strings.Fixed."*";
use Ada.Text_Io;
begin
for N in 1 .. Observations loop Histogram_20.Increment_Bucket (Modified_Random.Random); end loop;
Put ("Range Observations: "); Put (Observations'Image); Put (" Buckets: "); Put (Buckets'Image); New_Line; for I in Histogram_20.Bucket_Index'Range loop Put (Histogram_20.Range_Image (I)); Put (" "); Put (Column (Histogram_20.Observations_In (I) / Divider, Char)); New_Line; end loop;
end Modified_Distribution;</lang>
- Output:
Range Observations: 20000 Buckets: 25 0.00..0.04 **************************************************************************** 0.04..0.08 ************************************************************************ 0.08..0.12 ************************************************************** 0.12..0.16 ********************************************************** 0.16..0.20 ************************************************** 0.20..0.24 ****************************************** 0.24..0.28 *************************************** 0.28..0.32 ******************************* 0.32..0.36 ************************* 0.36..0.40 ******************* 0.40..0.44 ************ 0.44..0.48 ***** 0.48..0.52 * 0.52..0.56 ******* 0.56..0.60 *********** 0.60..0.64 ******************** 0.64..0.68 *************************** 0.68..0.72 ********************************** 0.72..0.76 ************************************** 0.76..0.80 ******************************************** 0.80..0.84 ************************************************* 0.84..0.88 ********************************************************** 0.88..0.92 **************************************************************** 0.92..0.96 ******************************************************************* 0.96..1.00 ************************************************************************
Factor
<lang factor>USING: assocs assocs.extras formatting io kernel math math.functions math.statistics random sequences tools.memory.private ;
- modifier ( x -- y ) 0.5 over 0.5 < [ swap ] when - dup + ;
- random-unit-by ( quot: ( x -- y ) -- z )
random-unit dup pick call random-unit 2dup > [ 2drop nip ] [ 3drop random-unit-by ] if ; inline recursive
- data ( n quot bins -- histogram )
'[ _ random-unit-by _ * >integer ] replicate histogram ; inline
- .histogram ( h -- )
h assoc-size :> buckets ! number of buckets h sum-values :> total ! items in histogram h values supremum :> max ! largest bucket (as in most occurrences) 40 :> size ! max size of a bar
total commas buckets "Bin Histogram (%s items, %d buckets)\n" printf
h [| k v | k buckets / dup buckets recip + "[%.2f, %.2f) " printf size v * max / ceiling [ "▇" write ] times bl bl v commas print ] assoc-each ;
"Modified random distribution of values in [0, 1):" print nl 100,000 [ modifier ] 20 data .histogram</lang>
- Output:
Modified random distribution of values in [0, 1): Bin Histogram (100,000 items, 20 buckets) [0.00, 0.05) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 9,416 [0.05, 0.10) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 8,498 [0.10, 0.15) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 7,432 [0.15, 0.20) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 6,415 [0.20, 0.25) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 5,558 [0.25, 0.30) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 4,489 [0.30, 0.35) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 3,538 [0.35, 0.40) ▇▇▇▇▇▇▇▇▇▇▇ 2,532 [0.40, 0.45) ▇▇▇▇▇▇▇ 1,553 [0.45, 0.50) ▇▇▇ 490 [0.50, 0.55) ▇▇▇ 517 [0.55, 0.60) ▇▇▇▇▇▇▇ 1,467 [0.60, 0.65) ▇▇▇▇▇▇▇▇▇▇▇ 2,519 [0.65, 0.70) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 3,559 [0.70, 0.75) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 4,546 [0.75, 0.80) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 5,569 [0.80, 0.85) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 6,444 [0.85, 0.90) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 7,428 [0.90, 0.95) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 8,487 [0.95, 1.00) ▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 9,543
FreeBASIC
<lang freebasic>#define NRUNS 100000
- define NBINS 20
function modifier( x as double ) as double
return iif(x < 0.5, 2*(0.5 - x), 2*(x - 0.5))
end function
function modrand() as double
dim as double random1, random2 do random1 = rnd random2 = rnd if random2 < modifier(random1) then return random1 endif loop
end function
function histo( byval bn as uinteger ) as string
dim as double db = NRUNS/(50*NBINS) dim as string h while bn > db: h = h + "#" bn -= db wend return h
end function
dim as uinteger bins(0 to NBINS-1), i, b dim as double db = 1./NBINS, rand
randomize timer
for i = 1 to NRUNS
rand = modrand() b = int(rand/db) bins(b) += 1
next i
for b = 0 to NBINS-1
print using "Bin ## (#.## to #.##): & ####";b;b*db;(b+1)*db;histo(bins(b));bins(b)
next b</lang>
- Output:
Bin 0 (0.00 to 0.05): ############################################################################################## 9479 Bin 1 (0.05 to 0.10): #################################################################################### 8499 Bin 2 (0.10 to 0.15): ########################################################################## 7416 Bin 3 (0.15 to 0.20): ################################################################## 6650 Bin 4 (0.20 to 0.25): ###################################################### 5457 Bin 5 (0.25 to 0.30): ############################################ 4416 Bin 6 (0.30 to 0.35): ################################## 3469 Bin 7 (0.35 to 0.40): ######################## 2481 Bin 8 (0.40 to 0.45): ############## 1466 Bin 9 (0.45 to 0.50): #### 489 Bin 10 (0.50 to 0.55): #### 475 Bin 11 (0.55 to 0.60): ############## 1472 Bin 12 (0.60 to 0.65): ######################### 2548 Bin 13 (0.65 to 0.70): #################################### 3617 Bin 14 (0.70 to 0.75): ############################################# 4538 Bin 15 (0.75 to 0.80): ####################################################### 5590 Bin 16 (0.80 to 0.85): ############################################################### 6395 Bin 17 (0.85 to 0.90): ########################################################################### 7538 Bin 18 (0.90 to 0.95): #################################################################################### 8401 Bin 19 (0.95 to 1.00): ################################################################################################ 9604
Go
<lang go>package main
import (
"fmt" "math" "math/rand" "strings" "time"
)
func rng(modifier func(x float64) float64) float64 {
for { r1 := rand.Float64() r2 := rand.Float64() if r2 < modifier(r1) { return r1 } }
}
func commatize(n int) string {
s := fmt.Sprintf("%d", n) if n < 0 { s = s[1:] } le := len(s) for i := le - 3; i >= 1; i -= 3 { s = s[0:i] + "," + s[i:] } if n >= 0 { return s } return "-" + s
}
func main() {
rand.Seed(time.Now().UnixNano()) modifier := func(x float64) float64 { if x < 0.5 { return 2 * (0.5 - x) } return 2 * (x - 0.5) } const ( N = 100000 NUM_BINS = 20 HIST_CHAR = "■" HIST_CHAR_SIZE = 125 ) bins := make([]int, NUM_BINS) // all zero by default binSize := 1.0 / NUM_BINS for i := 0; i < N; i++ { rn := rng(modifier) bn := int(math.Floor(rn / binSize)) bins[bn]++ }
fmt.Println("Modified random distribution with", commatize(N), "samples in range [0, 1):\n") fmt.Println(" Range Number of samples within that range") for i := 0; i < NUM_BINS; i++ { hist := strings.Repeat(HIST_CHAR, int(math.Round(float64(bins[i])/HIST_CHAR_SIZE))) fi := float64(i) fmt.Printf("%4.2f ..< %4.2f %s %s\n", binSize*fi, binSize*(fi+1), hist, commatize(bins[i])) }
}</lang>
- Output:
Specimen run:
Modified random distribution with 100,000 samples in range [0, 1): Range Number of samples within that range 0.00 ..< 0.05 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 9,396 0.05 ..< 0.10 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 8,434 0.10 ..< 0.15 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 7,484 0.15 ..< 0.20 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 6,576 0.20 ..< 0.25 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 5,516 0.25 ..< 0.30 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 4,625 0.30 ..< 0.35 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 3,478 0.35 ..< 0.40 ■■■■■■■■■■■■■■■■■■■■ 2,506 0.40 ..< 0.45 ■■■■■■■■■■■■ 1,504 0.45 ..< 0.50 ■■■■ 505 0.50 ..< 0.55 ■■■■ 511 0.55 ..< 0.60 ■■■■■■■■■■■■■ 1,563 0.60 ..< 0.65 ■■■■■■■■■■■■■■■■■■■■■ 2,582 0.65 ..< 0.70 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 3,520 0.70 ..< 0.75 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 4,326 0.75 ..< 0.80 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 5,489 0.80 ..< 0.85 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 6,589 0.85 ..< 0.90 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 7,472 0.90 ..< 0.95 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 8,592 0.95 ..< 1.00 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 9,332
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
Mathematica/Wolfram Language
<lang Mathematica>ClearAll[Modifier, CreateRandomNumber] Modifier[x_] := If[x < 0.5, 2 (0.5 - x), 2 (x - 0.5)] CreateRandomNumber[] := Module[{r1, r2, done = True},
While[done, r1 = RandomReal[]; r2 = RandomReal[]; If[r2 < Modifier[r1], Return[r1]; done = False ] ] ]
numbers = Table[CreateRandomNumber[], 100000]; {bins, counts} = HistogramList[numbers, {0, 1, 0.05}, "PDF"]; Grid[MapThread[{#1, " - ", StringJoin@ConstantArray["X", Round[20 #2]]} &, {Partition[bins, 2, 1], counts}], Alignment -> Left]</lang>
- Output:
{0.,0.05} - XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX {0.05,0.1} - XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX {0.1,0.15} - XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX {0.15,0.2} - XXXXXXXXXXXXXXXXXXXXXXXXXX {0.2,0.25} - XXXXXXXXXXXXXXXXXXXXXX {0.25,0.3} - XXXXXXXXXXXXXXXXXX {0.3,0.35} - XXXXXXXXXXXXXX {0.35,0.4} - XXXXXXXXXX {0.4,0.45} - XXXXXX {0.45,0.5} - XX {0.5,0.55} - XX {0.55,0.6} - XXXXXX {0.6,0.65} - XXXXXXXXXX {0.65,0.7} - XXXXXXXXXXXXXX {0.7,0.75} - XXXXXXXXXXXXXXXXXX {0.75,0.8} - XXXXXXXXXXXXXXXXXXXXXX {0.8,0.85} - XXXXXXXXXXXXXXXXXXXXXXXXXX {0.85,0.9} - XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX {0.9,0.95} - XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX {0.95,1.} - XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
Nim
<lang Nim>import random, strformat, strutils, sugar
type ValRange = range[0.0..1.0]
func modifier(x: ValRange): ValRange =
if x < 0.5: 2 * (0.5 - x) else: 2 * (x - 0.5)
proc rand(modifier: (float) -> float): ValRange =
while true: let r1 = rand(1.0) let r2 = rand(1.0) if r2 < modifier(r1): return r1
const
N = 100_000 NumBins = 20 HistChar = "■" HistCharSize = 125 BinSize = 1 / NumBins
randomize()
var bins: array[NumBins, int] for i in 0..<N:
let rn = rand(modifier) let bn = int(rn / BinSize) inc bins[bn]
echo &"Modified random distribution with {N} samples in range [0, 1):" echo " Range Number of samples within that range" for i in 0..<NumBins:
let hist = repeat(HistChar, (bins[i] / HistCharSize).toInt) echo &"{BinSize * float(i):4.2f} ..< {BinSize * float(i + 1):4.2f} {hist} {bins[i]}"</lang>
- Output:
Modified random distribution with 100000 samples in range [0, 1): Range Number of samples within that range 0.00 ..< 0.05 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 9480 0.05 ..< 0.10 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 8469 0.10 ..< 0.15 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 7631 0.15 ..< 0.20 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 6484 0.20 ..< 0.25 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 5472 0.25 ..< 0.30 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 4327 0.30 ..< 0.35 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 3523 0.35 ..< 0.40 ■■■■■■■■■■■■■■■■■■■■ 2528 0.40 ..< 0.45 ■■■■■■■■■■■■ 1500 0.45 ..< 0.50 ■■■■ 444 0.50 ..< 0.55 ■■■■ 513 0.55 ..< 0.60 ■■■■■■■■■■■■ 1536 0.60 ..< 0.65 ■■■■■■■■■■■■■■■■■■■■ 2459 0.65 ..< 0.70 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 3505 0.70 ..< 0.75 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 4600 0.75 ..< 0.80 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 5525 0.80 ..< 0.85 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 6512 0.85 ..< 0.90 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 7482 0.90 ..< 0.95 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 8581 0.95 ..< 1.00 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 9429
Perl
Uses any supplied distribution function, but defaults to uniform otherwise. <lang perl>use strict; use warnings; use List::Util 'max';
sub distribution {
my %param = ( function => \&{scalar sub {return 1}}, sample_size => 1e5, @_); my @values; do { my($r1, $r2) = (rand, rand); push @values, $r1 if &{$param{function}}($r1) > $r2; } until @values == $param{sample_size}; wantarray ? @values : \@values;
}
sub modifier_notch {
my($x) = @_; return 2 * ( $x < 1/2 ? ( 1/2 - $x ) : ( $x - 1/2 ) );
}
sub print_histogram {
our %param = (n_bins => 10, width => 80, @_); my %counts; $counts{ int($_ * $param{n_bins}) / $param{n_bins} }++ for @{$param{data}}; our $max_value = max values %counts; print "Bin Counts Histogram\n"; printf "%4.2f %6d: %s\n", $_, $counts{$_}, hist($counts{$_}) for sort keys %counts; sub hist { scalar ('■') x ( $param{width} * $_[0] / $max_value ) }
}
print_histogram( data => \@{ distribution() } ); print "\n\n";
my @samples = distribution( function => \&modifier_notch, sample_size => 50_000); print_histogram( data => \@samples, n_bins => 20, width => 64);</lang>
- Output:
Bin Counts Histogram 0.00 10114: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.10 9958: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.20 9960: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.30 10043: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.40 9874: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.50 10013: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.60 10085: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.70 9877: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.80 10079: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.90 9997: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ Bin Counts Histogram 0.00 4772: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.05 4329: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.10 3728: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.15 3249: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.20 2749: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.25 2163: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.30 1735: ■■■■■■■■■■■■■■■■■■■■■■■ 0.35 1317: ■■■■■■■■■■■■■■■■■ 0.40 764: ■■■■■■■■■■ 0.45 259: ■■■ 0.50 231: ■■■ 0.55 721: ■■■■■■■■■ 0.60 1255: ■■■■■■■■■■■■■■■■ 0.65 1730: ■■■■■■■■■■■■■■■■■■■■■■■ 0.70 2282: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.75 2720: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.80 3302: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.85 3712: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.90 4219: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.95 4763: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
Phix
<lang Phix>function rng(integer modifier)
while true do atom r1 := rnd() if rnd() < modifier(r1) then return r1 end if end while
end function
function modifier(atom x)
return iff(x < 0.5 ? 2 * (0.5 - x) : 2 * (x - 0.5))
end function
constant N = 100000,
NUM_BINS = 20, HIST_CHAR_SIZE = 125, BIN_SIZE = 1/NUM_BINS, LO = sq_mul(tagset(NUM_BINS-1,0),BIN_SIZE), HI = sq_mul(tagset(NUM_BINS),BIN_SIZE), LBLS = apply(true,sprintf,{{"[%4.2f,%4.2f)"},columnize({LO,HI})})
sequence bins := repeat(0, NUM_BINS) for i=1 to N do
bins[floor(rng(modifier) / BIN_SIZE)+1] += 1
end for
printf(1,"Modified random distribution with %,d samples in range [0, 1):\n\n",N) for i=1 to NUM_BINS do
sequence hist := repeat('#', round(bins[i]/HIST_CHAR_SIZE)) printf(1,"%s %s %,d\n", {LBLS[i], hist, bins[i]})
end for</lang>
- Output:
Modified random distribution with 100,000 samples in range [0, 1): [0.00,0.05) ############################################################################ 9,521 [0.05,0.10) #################################################################### 8,449 [0.10,0.15) ############################################################ 7,519 [0.15,0.20) ##################################################### 6,651 [0.20,0.25) ############################################ 5,470 [0.25,0.30) #################################### 4,504 [0.30,0.35) ########################### 3,364 [0.35,0.40) #################### 2,475 [0.40,0.45) ############ 1,494 [0.45,0.50) #### 518 [0.50,0.55) #### 482 [0.55,0.60) ############ 1,536 [0.60,0.65) ##################### 2,568 [0.65,0.70) ############################ 3,498 [0.70,0.75) #################################### 4,559 [0.75,0.80) ############################################ 5,447 [0.80,0.85) #################################################### 6,512 [0.85,0.90) ############################################################ 7,486 [0.90,0.95) #################################################################### 8,484 [0.95,1.00) ############################################################################ 9,463
plot
A simple graphical plot. Note the labels are on the X-axis, so it's v
-shaped, not <
-shaped: IupPlot does not support putting user-supplied labels on the Y-axis.
<lang Phix>include pGUI.e
IupOpen()
Ihandle plot = IupPlot("GRID=YES, AXS_YAUTOMIN=NO")
IupPlotBegin(plot, true) -- (true means x-axis are labels)
for i=1 to length(bins) do
IupPlotAddStr(plot, LBLS[i], bins[i]);
end for {} = IupPlotEnd(plot) IupSetAttribute(plot,"DS_MODE","BAR") IupSetAttribute(plot,"DS_COLOR",IUP_DARK_BLUE) IupShow(IupDialog(plot, `TITLE=Histogram, RASTERSIZE=1300x850`)) IupMainLoop() IupClose()</lang>
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
If a vertical histograph (instead of a < shaped horizontal histograph) were to be used, it would be a V shaped. <lang rexx>/*REXX program generates a "<" 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
Raku
<lang perl6>sub modified_random_distribution ( Code $modifier --> Seq ) {
return lazy gather loop { my ( $r1, $r2 ) = rand, rand; take $r1 if $modifier.($r1) > $r2; }
} sub modifier ( Numeric $x --> Numeric ) {
return 2 * ( $x < 1/2 ?? ( 1/2 - $x ) !! ( $x - 1/2 ) );
} sub print_histogram ( @data, :$n-bins, :$width ) { # Assumes minimum of zero.
my %counts = bag @data.map: { floor( $_ * $n-bins ) / $n-bins }; my $max_value = %counts.values.max; sub hist { '■' x ( $width * $^count / $max_value ) } say ' Bin, Counts: Histogram'; printf "%4.2f, %6d: %s\n", .key, .value, hist(.value) for %counts.sort;
}
my @d = modified_random_distribution( &modifier );
print_histogram( @d.head(50_000), :n-bins(20), :width(64) );</lang>
- Output:
Bin, Counts: Histogram 0.00, 4718: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.05, 4346: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.10, 3685: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.15, 3246: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.20, 2734: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.25, 2359: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.30, 1702: ■■■■■■■■■■■■■■■■■■■■■■ 0.35, 1283: ■■■■■■■■■■■■■■■■■ 0.40, 702: ■■■■■■■■■ 0.45, 250: ■■■ 0.50, 273: ■■■ 0.55, 745: ■■■■■■■■■■ 0.60, 1231: ■■■■■■■■■■■■■■■■ 0.65, 1757: ■■■■■■■■■■■■■■■■■■■■■■■ 0.70, 2209: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.75, 2738: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.80, 3255: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.85, 3741: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.90, 4268: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.95, 4758: ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■
Wren
<lang ecmascript>import "random" for Random import "/fmt" for Fmt
var rgen = Random.new()
var rng = Fn.new { |modifier|
while (true) { var r1 = rgen.float() var r2 = rgen.float() if (r2 < modifier.call(r1)) { return r1 } }
}
var modifier = Fn.new { |x| (x < 0.5) ? 2 * (0.5 - x) : 2 * (x - 0.5) }
var N = 100000 var NUM_BINS = 20 var HIST_CHAR = "■" var HIST_CHAR_SIZE = 125 var bins = List.filled(NUM_BINS, 0) var binSize = 1 / NUM_BINS for (i in 0...N) {
var rn = rng.call(modifier) var bn = (rn / binSize).floor bins[bn] = bins[bn] + 1
}
Fmt.print("Modified random distribution with $,d samples in range [0, 1):\n", N) System.print(" Range Number of samples within that range") for (i in 0...NUM_BINS) {
var hist = HIST_CHAR * (bins[i] / HIST_CHAR_SIZE).round Fmt.print("$4.2f ..< $4.2f $s $,d", binSize * i, binSize * (i + 1), hist, bins[i])
}</lang>
- Output:
Specimen run:
Modified random distribution with 100,000 samples in range [0, 1): Range Number of samples within that range 0.00 ..< 0.05 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 9,605 0.05 ..< 0.10 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 8,573 0.10 ..< 0.15 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 7,440 0.15 ..< 0.20 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 6,582 0.20 ..< 0.25 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 5,482 0.25 ..< 0.30 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 4,472 0.30 ..< 0.35 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 3,478 0.35 ..< 0.40 ■■■■■■■■■■■■■■■■■■■■ 2,497 0.40 ..< 0.45 ■■■■■■■■■■■■ 1,519 0.45 ..< 0.50 ■■■■ 489 0.50 ..< 0.55 ■■■■ 485 0.55 ..< 0.60 ■■■■■■■■■■■■ 1,453 0.60 ..< 0.65 ■■■■■■■■■■■■■■■■■■■■ 2,477 0.65 ..< 0.70 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 3,492 0.70 ..< 0.75 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 4,453 0.75 ..< 0.80 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 5,535 0.80 ..< 0.85 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 6,480 0.85 ..< 0.90 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 7,573 0.90 ..< 0.95 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 8,372 0.95 ..< 1.00 ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 9,543