Unbias a random generator: Difference between revisions
(→{{header|Haskell}}: Changed to universal monadic solution) |
(→{{header|Haskell}}: Added examples) |
||
| Line 681: | Line 681: | ||
=={{header|Haskell}}== |
=={{header|Haskell}}== |
||
The first task: |
|||
<lang haskell>import Control.Monad.Random |
<lang haskell>import Control.Monad.Random |
||
import Control.Monad |
import Control.Monad |
||
| Line 687: | Line 687: | ||
randN :: MonadRandom m => Int -> m Int |
randN :: MonadRandom m => Int -> m Int |
||
randN n = fromList [(0, fromIntegral n-1), (1, 1)] |
randN n = fromList [(0, fromIntegral n-1), (1, 1)]</lang> |
||
Examples of use: |
|||
-- works for any random generator |
|||
<pre>λ> replicateM 20 (randN 2) |
|||
| ⚫ | |||
[0,0,1,0,0,1,0,1,1,0,0,1,1,1,1,0,1,1,0,0] |
|||
| ⚫ | |||
λ> replicateM 20 (randN 5) |
|||
[0,1,0,0,0,0,0,0,1,0,0,0,0,1,0,0,1,0,1,0]</pre> |
|||
The second task. Returns the unbiased generator for any given random generator. |
|||
| ⚫ | |||
| ⚫ | |||
y <- g |
y <- g |
||
if x /= y then return y else unbias g |
if x /= y then return y else unbias g</lang> |
||
Examples of use: |
|||
<pre>λ> replicateM 20 (unbiased (randN 5)) |
|||
[0,0,1,0,1,1,1,0,0,0,1,1,1,0,1,1,0,0,1,0] |
|||
λ> replicateM 20 (unbias (fromList [(True,10),(False,1)])) |
|||
[True,True,False,True,True,True,False,True,False,True,True,False,False,True,False,True,True,False,False,True]</pre> |
|||
The third task: |
|||
main = forM_ [3..6] showCounts |
<lang Haskell>main = forM_ [3..6] showCounts |
||
where |
where |
||
showCounts b = do |
showCounts b = do |
||
Revision as of 12:07, 23 October 2016
You are encouraged to solve this task according to the task description, using any language you may know.
Given a weighted one bit generator of random numbers where the probability of a one occuring, , is not the same as , the probability of a zero occuring, the probability of the occurrence of a one followed by a zero is × . This is the same as the probability of a zero followed by a one: × .
- Task details
- Use your language's random number generator to create a function/method/subroutine/... randN that returns a one or a zero, but with one occurring, on average, 1 out of N times, where N is an integer from the range 3 to 6 inclusive.
- Create a function unbiased that uses only randN as its source of randomness to become an unbiased generator of random ones and zeroes.
- For N over its range, generate and show counts of the outputs of randN and unbiased(randN).
The actual unbiasing should be done by generating two numbers at a time from randN and only returning a 1 or 0 if they are different. As long as you always return the first number or always return the second number, the probabilities discussed above should take over the biased probability of randN.
This task is an implementation of Von Neumann debiasing, first described in a 1951 paper.
Ada
<lang Ada>with Ada.Text_IO; with Ada.Numerics.Discrete_Random;
procedure Bias_Unbias is
Modulus: constant Integer := 60; -- lcm of {3,4,5,6}
type M is mod Modulus;
package Rand is new Ada.Numerics.Discrete_Random(M);
Gen: Rand.Generator;
subtype Bit is Integer range 0 .. 1;
function Biased_Bit(Bias_Base: Integer) return Bit is
begin
if (Integer(Rand.Random(Gen))* Bias_Base) / Modulus > 0 then
return 0;
else
return 1;
end if;
end Biased_Bit;
function Unbiased_Bit(Bias_Base: Integer) return Bit is
A, B: Bit := 0;
begin
while A = B loop
A := Biased_Bit(Bias_Base);
B := Biased_Bit(Bias_Base);
end loop;
return A;
end Unbiased_Bit;
package FIO is new Ada.Text_IO.Float_IO(Float);
Counter_B, Counter_U: Natural; Number_Of_Samples: constant Natural := 10_000;
begin
Rand.Reset(Gen);
Ada.Text_IO.Put_Line(" I Biased% UnBiased%");
for I in 3 .. 6 loop
Counter_B := 0;
Counter_U := 0;
for J in 1 .. Number_Of_Samples loop
Counter_B := Counter_B + Biased_Bit(I);
Counter_U := Counter_U + Unbiased_Bit(I);
end loop;
Ada.Text_IO.Put(Integer'Image(I));
FIO.Put(100.0 * Float(Counter_B) / Float(Number_Of_Samples), 5, 2, 0);
FIO.Put(100.0 * Float(Counter_U) / Float(Number_Of_Samples), 5, 2, 0);
Ada.Text_IO.New_Line;
end loop;
end Bias_Unbias;</lang>
Output:
I Biased% UnBiased% 3 32.87 49.80 4 24.49 50.22 5 19.73 50.05 6 16.75 50.19
Aime
<lang aime>integer biased(integer bias) {
return 1 ^ min(drand(bias - 1), 1);
}
integer unbiased(integer bias) {
integer a;
while ((a = biased(bias)) == biased(bias)) {
}
return a;
}
integer main(void) {
integer b, n, cb, cu, i;
n = 10000;
b = 3;
while (b <= 6) {
i = 0;
cb = 0;
cu = 0;
while (i < n) {
cb += biased(b);
cu += unbiased(b);
i += 1;
}
o_form("bias ~: /d2p2/%% vs /d2p2/%%\n", b, 100r * cb / n,
100r * cu / n);
b += 1; }
return 0;
}</lang>
Output:
bias 3: 33.51% vs 50.27% bias 4: 24.97% vs 49.99% bias 5: 19.93% vs 49.92% bias 6: 16.32% vs 49.36%
AutoHotkey
<lang AutoHotkey>Biased(){
Random, q, 0, 4 return q=4
} Unbiased(){
Loop
If ((a := Biased()) != biased())
return a
} Loop 1000
t .= biased(), t2 .= unbiased()
StringReplace, junk, t2, 1, , UseErrorLevel MsgBox % "Unbiased probability of a 1 occurring: " Errorlevel/1000 StringReplace, junk, t, 1, , UseErrorLevel MsgBox % "biased probability of a 1 occurring: " Errorlevel/1000</lang>
BBC BASIC
<lang bbcbasic> FOR N% = 3 TO 6
biased% = 0
unbiased% = 0
FOR I% = 1 TO 10000
IF FNrandN(N%) biased% += 1
IF FNunbiased(N%) unbiased% += 1
NEXT
PRINT "N = ";N% " : biased = "; biased%/100 "%, unbiased = "; unbiased%/100 "%"
NEXT
END
DEF FNunbiased(N%)
LOCAL A%,B%
REPEAT
A% = FNrandN(N%)
B% = FNrandN(N%)
UNTIL A%<>B%
= A%
DEF FNrandN(N%) = -(RND(N%) = 1)</lang>
Output:
N = 3 : biased = 33.57%, unbiased = 49.94% N = 4 : biased = 25.34%, unbiased = 50.76% N = 5 : biased = 20.06%, unbiased = 50.04% N = 6 : biased = 16.25%, unbiased = 50.13%
C
<lang C>#include <stdio.h>
- include <stdlib.h>
int biased(int bias) { /* balance out the bins, being pedantic */ int r, rand_max = RAND_MAX - (RAND_MAX % bias); while ((r = rand()) > rand_max); return r < rand_max / bias; }
int unbiased(int bias) { int a; while ((a = biased(bias)) == biased(bias)); return a; }
int main() { int b, n = 10000, cb, cu, i; for (b = 3; b <= 6; b++) { for (i = cb = cu = 0; i < n; i++) { cb += biased(b); cu += unbiased(b); } printf("bias %d: %5.3f%% vs %5.3f%%\n", b, 100. * cb / n, 100. * cu / n); }
return 0; }</lang> output
bias 3: 33.090% vs 49.710% bias 4: 25.130% vs 49.430% bias 5: 19.760% vs 49.650% bias 6: 16.740% vs 50.030%
C#
<lang c sharp>using System;
namespace Unbias {
internal class Program
{
private static void Main(string[] args)
{
// Demonstrate.
for (int n = 3; n <= 6; n++)
{
int biasedZero = 0, biasedOne = 0, unbiasedZero = 0, unbiasedOne = 0;
for (int i = 0; i < 100000; i++)
{
if (randN(n))
biasedOne++;
else
biasedZero++;
if (Unbiased(n))
unbiasedOne++;
else
unbiasedZero++;
}
Console.WriteLine("(N = {0}):".PadRight(17) + "# of 0\t# of 1\t% of 0\t% of 1", n);
Console.WriteLine("Biased:".PadRight(15) + "{0}\t{1}\t{2}\t{3}",
biasedZero, biasedOne,
biasedZero/1000, biasedOne/1000);
Console.WriteLine("Unbiased:".PadRight(15) + "{0}\t{1}\t{2}\t{3}",
unbiasedZero, unbiasedOne,
unbiasedZero/1000, unbiasedOne/1000);
}
}
private static bool Unbiased(int n)
{
bool flip1, flip2;
/* Flip twice, and check if the values are the same.
* If so, flip again. Otherwise, return the value of the first flip. */
do
{
flip1 = randN(n);
flip2 = randN(n);
} while (flip1 == flip2);
return flip1;
}
private static readonly Random random = new Random();
private static bool randN(int n)
{
// Has an 1/n chance of returning 1. Otherwise it returns 0.
return random.Next(0, n) == 0;
}
}
}</lang>
Sample Output
(N = 3): # of 0 # of 1 % of 0 % of 1 Biased: 66867 33133 66 33 Unbiased: 49843 50157 49 50 (N = 4): # of 0 # of 1 % of 0 % of 1 Biased: 74942 25058 74 25 Unbiased: 50192 49808 50 49 (N = 5): # of 0 # of 1 % of 0 % of 1 Biased: 80203 19797 80 19 Unbiased: 49928 50072 49 50 (N = 6): # of 0 # of 1 % of 0 % of 1 Biased: 83205 16795 83 16 Unbiased: 49744 50256 49 50
Clojure
<lang Clojure>(defn biased [n]
(if (< (rand 2) (/ n)) 0 1))
(defn unbiased [n]
(loop [a 0 b 0]
(if (= a b)
(recur (biased n) (biased n))
a)))
(for [n (range 3 7)]
[n (double (/ (apply + (take 50000 (repeatedly #(biased n)))) 50000)) (double (/ (apply + (take 50000 (repeatedly #(unbiased n)))) 50000))])
([3 0.83292 0.50422]
[4 0.87684 0.5023] [5 0.90122 0.49728] [6 0.91526 0.5])</lang>
CoffeeScript
<lang coffeescript> biased_rand_function = (n) ->
# return a function that returns 0/1 with
# 1 appearing only 1/Nth of the time
cap = 1/n
->
if Math.random() < cap
1
else
0
unbiased_function = (f) ->
->
while true
[n1, n2] = [f(), f()]
return n1 if n1 + n2 == 1
stats = (label, f) ->
cnt = 0
sample_size = 10000000
for i in [1...sample_size]
cnt += 1 if f() == 1
console.log "ratio of 1s: #{cnt / sample_size} [#{label}]"
for n in [3..6]
console.log "\n---------- n = #{n}"
f_biased = biased_rand_function(n)
f_unbiased = unbiased_function f_biased
stats "biased", f_biased
stats "unbiased", f_unbiased
</lang> output
> coffee unbiased.coffee ---------- n = 3 ratio of 1s: 0.3333343 [biased] ratio of 1s: 0.4999514 [unbiased] ---------- n = 4 ratio of 1s: 0.2499751 [biased] ratio of 1s: 0.4998067 [unbiased] ---------- n = 5 ratio of 1s: 0.199729 [biased] ratio of 1s: 0.5003183 [unbiased] ---------- n = 6 ratio of 1s: 0.1664843 [biased] ratio of 1s: 0.4997813 [unbiased]
Common Lisp
<lang lisp>(defun biased (n) (if (zerop (random n)) 0 1))
(defun unbiased (n)
(loop with x do
(if (/= (setf x (biased n)) (biased n))
(return x))))
(loop for n from 3 to 6 do
(let ((u (loop repeat 10000 collect (unbiased n)))
(b (loop repeat 10000 collect (biased n)))) (format t "~a: unbiased ~d biased ~d~%" n (count 0 u) (count 0 b))))</lang> output
3: unbiased 4992 biased 3361 4: unbiased 4988 biased 2472 5: unbiased 5019 biased 1987 6: unbiased 4913 biased 1658
D
<lang d>import std.stdio, std.random, std.algorithm, std.range, std.functional;
enum biased = (in int n) /*nothrow*/ => uniform01 < (1.0 / n);
int unbiased(in int bias) /*nothrow*/ {
int a;
while ((a = bias.biased) == bias.biased) {}
return a;
}
void main() {
enum M = 500_000;
foreach (immutable n; 3 .. 7)
writefln("%d: %2.3f%% %2.3f%%", n,
M.iota.map!(_=> n.biased).sum * 100.0 / M,
M.iota.map!(_=> n.unbiased).sum * 100.0 / M);
}</lang>
- Output:
3: 33.441% 49.964% 4: 24.953% 49.910% 5: 19.958% 49.987% 6: 16.660% 49.890%
Elixir
<lang elixir>defmodule Random do
def randN(n) do
if :rand.uniform(n) == 1, do: 1, else: 0
end
def unbiased(n) do
{x, y} = {randN(n), randN(n)}
if x != y, do: x, else: unbiased(n)
end
end
IO.puts "N biased unbiased" m = 10000 for n <- 3..6 do
xs = for _ <- 1..m, do: Random.randN(n)
ys = for _ <- 1..m, do: Random.unbiased(n)
IO.puts "#{n} #{Enum.sum(xs) / m} #{Enum.sum(ys) / m}"
end</lang>
- Output:
N biased unbiased 3 0.3356 0.5043 4 0.2523 0.4996 5 0.2027 0.5041 6 0.1647 0.4912
ERRE
<lang ERRE>PROGRAM UNBIAS
FUNCTION RANDN(N)
RANDN=INT(1+N*RND(1))=1
END FUNCTION
PROCEDURE UNBIASED(N->RIS)
LOCAL A,B
REPEAT
A=RANDN(N)
B=RANDN(N)
UNTIL A<>B
RIS=A
END PROCEDURE
BEGIN
PRINT(CHR$(12);) ! CLS RANDOMIZE(TIMER)
FOR N=3 TO 6 DO
BIASED=0
UNBIASED=0
FOR I=1 TO 10000 DO
IF RANDN(N) THEN biased+=1
UNBIASED(N->RIS)
IF RIS THEN unbiased+=+1
END FOR
PRINT("N =";N;" : biased =";biased/100;", unbiased =";unbiased/100)
END FOR
END PROGRAM </lang>
- Output:
N = 3 : biased = 32.66 , unbiased = 49.14 N = 4 : biased = 25.49 , unbiased = 49.92 N = 5 : biased = 20.53 , unbiased = 50 N = 6 : biased = 17.43 , unbiased = 50.43
Euphoria
<lang euphoria>function randN(integer N)
return rand(N) = 1
end function
function unbiased(integer N)
integer a
while 1 do
a = randN(N)
if a != randN(N) then
return a
end if
end while
end function
constant n = 10000 integer cb, cu for b = 3 to 6 do
cb = 0
cu = 0
for i = 1 to n do
cb += randN(b)
cu += unbiased(b)
end for
printf(1, "%d: %5.2f%% %5.2f%%\n", {b, 100 * cb / n, 100 * cu / n})
end for</lang>
Output:
3: 33.68% 49.94% 4: 24.93% 50.48% 5: 20.32% 49.97% 6: 16.98% 50.05%
F#
<lang fsharp>open System
let random = Random()
let randN = random.Next >> (=)0 >> Convert.ToInt32
let rec unbiased n =
let a = randN n if a <> randN n then a else unbiased n
[<EntryPoint>] let main argv =
let n = if argv.Length > 0 then UInt32.Parse(argv.[0]) |> int else 100000
for b = 3 to 6 do
let cb = ref 0
let cu = ref 0
for i = 1 to n do
cb := !cb + randN b
cu := !cu + unbiased b
printfn "%d: %5.2f%% %5.2f%%"
b (100. * float !cb / float n) (100. * float !cu / float n)
0</lang>
- Output:
3: 33.26% 49.97% 4: 25.02% 50.22% 5: 19.98% 50.00% 6: 16.64% 49.69%
Fortran
<lang fortran>program Bias_Unbias
implicit none
integer, parameter :: samples = 1000000
integer :: i, j
integer :: c1, c2, rand
do i = 3, 6
c1 = 0
c2 = 0
do j = 1, samples
rand = bias(i)
if (rand == 1) c1 = c1 + 1
rand = unbias(i)
if (rand == 1) c2 = c2 + 1
end do
write(*, "(i2,a,f8.3,a,f8.3,a)") i, ":", real(c1) * 100.0 / real(samples), &
"%", real(c2) * 100.0 / real(samples), "%"
end do
contains
function bias(n)
integer :: bias integer, intent(in) :: n real :: r
call random_number(r) if (r > 1 / real(n)) then bias = 0 else bias = 1 end if
end function
function unbias(n)
integer :: unbias integer, intent(in) :: n integer :: a, b
do a = bias(n) b = bias(n) if (a /= b) exit end do unbias = a
end function
end program</lang> Output:
3: 33.337% 49.971% 4: 24.945% 49.944% 5: 19.971% 49.987% 6: 16.688% 50.097%
GAP
<lang gap>RandNGen := function(n) local v, rand; v := [1 .. n - 1]*0; Add(v, 1); rand := function() return Random(v); end; return rand; end;
UnbiasedGen := function(rand) local unbiased; unbiased := function() local a, b; while true do a := rand(); b := rand(); if a <> b then break; fi; od; return a; end; return unbiased; end;
range := [2 .. 6]; v := List(range, RandNGen); w := List(v, UnbiasedGen); apply := gen -> Sum([1 .. 1000000], n -> gen());
- Some tests (2 is added as a witness, since in this case RandN is already unbiased)
PrintArray(TransposedMat([range, List(v, apply), List(w, apply)]));
- [ [ 2, 499991, 499041 ],
- [ 3, 333310, 500044 ],
- [ 4, 249851, 500663 ],
- [ 5, 200532, 500448 ],
- [ 6, 166746, 499859 ] ]</lang>
Go
<lang go>package main
import (
"fmt" "math/rand"
)
const samples = 1e6
func main() {
fmt.Println("Generator 1 count 0 count % 1 count")
for n := 3; n <= 6; n++ {
// function randN, per task description
randN := func() int {
if rand.Intn(n) == 0 {
return 1
}
return 0
}
var b [2]int
for x := 0; x < samples; x++ {
b[randN()]++
}
fmt.Printf("randN(%d) %7d %7d %5.2f%%\n",
n, b[1], b[0], float64(b[1])*100/samples)
// function unbiased, per task description
unbiased := func() (b int) {
for b = randN(); b == randN(); b = randN() {
}
return
}
var u [2]int
for x := 0; x < samples; x++ {
u[unbiased()]++
}
fmt.Printf("unbiased %7d %7d %5.2f%%\n",
u[1], u[0], float64(u[1])*100/samples)
}
}</lang> Output:
Generator 1 count 0 count % 1 count randN(3) 332711 667289 33.27% unbiased 499649 500351 49.96% randN(4) 249742 750258 24.97% unbiased 499434 500566 49.94% randN(5) 200318 799682 20.03% unbiased 499100 500900 49.91% randN(6) 166900 833100 16.69% unbiased 499973 500027 50.00%
Haskell
The first task: <lang haskell>import Control.Monad.Random import Control.Monad import Text.Printf
randN :: MonadRandom m => Int -> m Int randN n = fromList [(0, fromIntegral n-1), (1, 1)]</lang>
Examples of use:
λ> replicateM 20 (randN 2) [0,0,1,0,0,1,0,1,1,0,0,1,1,1,1,0,1,1,0,0] λ> replicateM 20 (randN 5) [0,1,0,0,0,0,0,0,1,0,0,0,0,1,0,0,1,0,1,0]
The second task. Returns the unbiased generator for any given random generator. <lang Haskell>unbiased :: (MonadRandom m, Eq x) => m x -> m x unbiased g = do x <- g
y <- g
if x /= y then return y else unbias g</lang>
Examples of use:
λ> replicateM 20 (unbiased (randN 5)) [0,0,1,0,1,1,1,0,0,0,1,1,1,0,1,1,0,0,1,0] λ> replicateM 20 (unbias (fromList [(True,10),(False,1)])) [True,True,False,True,True,True,False,True,False,True,True,False,False,True,False,True,True,False,False,True]
The third task: <lang Haskell>main = forM_ [3..6] showCounts
where
showCounts b = do
r1 <- counts (randN b)
r2 <- counts (unbias (randN b))
printf "n = %d biased: %d%% unbiased: %d%%\n" b r1 r2
counts g = (`div` 100) . length . filter (== 1) <$> replicateM 10000 g</lang>
Output:
n = 3 biased: 33% unbiased: 49% n = 4 biased: 24% unbiased: 50% n = 5 biased: 19% unbiased: 50% n = 6 biased: 16% unbiased: 49%
Icon and Unicon
This solution works in both languages. Both randN and unbiased are generators in the Icon/Unicon sense. <lang unicon>procedure main(A)
iters := \A[1] | 10000
write("ratios of 0 to 1 from ",iters," trials:")
every n := 3 to 6 do {
results_randN := table(0)
results_unbiased := table(0)
every 1 to iters do {
results_randN[randN(n)] +:= 1
results_unbiased[unbiased(n)] +:= 1
}
showResults(n, "randN", results_randN)
showResults(n, "unbiased", results_unbiased)
}
end
procedure showResults(n, s, t)
write(n," ",left(s,9),":",t[0],"/",t[1]," = ",t[0]/real(t[1]))
end
procedure unbiased(n)
repeat {
n1 := randN(n)
n2 := randN(n)
if n1 ~= n2 then suspend n1
}
end
procedure randN(n)
repeat suspend if 1 = ?n then 1 else 0
end</lang> and a sample run:
->ubrn 100000 ratios of 0 to 1 from 100000 trials: 3 randN :66804/33196 = 2.012411133871551 3 unbiased :49812/50188 = 0.9925081692834941 4 randN :75017/24983 = 3.002721850858584 4 unbiased :50000/50000 = 1.0 5 randN :79990/20010 = 3.997501249375312 5 unbiased :50073/49927 = 1.002924269433373 6 randN :83305/16695 = 4.989817310572027 6 unbiased :49911/50089 = 0.9964463255405378 ->
J
<lang j>randN=: 0 = ? unbiased=: i.@# { ::$: 2 | 0 3 -.~ _2 #.\ 4&* randN@# ]</lang>
Example use:
<lang j> randN 10#6 1 0 0 0 1 0 0 0 0 0
unbiased 10#6
1 0 0 1 0 0 1 0 1 1</lang>
Some example counts (these are counts of the number of 1s which appear in a test involving 100 random numbers):
<lang j> +/randN 100#3 30
+/randN 100#4
20
+/randN 100#5
18
+/randN 100#6
18
+/unbiased 100#3
49
+/unbiased 100#4
46
+/unbiased 100#5
49
+/unbiased 100#6
47</lang>
Note that these results are random. For example, a re-run of +/randN 100#5 gave 25 as its result, and a re-run of +/unbiased 100#5 gave 52 as its result.
Java
<lang java>public class Bias {
public static boolean biased(int n) {
return Math.random() < 1.0 / n;
}
public static boolean unbiased(int n) {
boolean a, b;
do {
a = biased(n);
b = biased(n);
} while (a == b);
return a;
}
public static void main(String[] args) {
final int M = 50000;
for (int n = 3; n < 7; n++) {
int c1 = 0, c2 = 0;
for (int i = 0; i < M; i++) {
c1 += biased(n) ? 1 : 0;
c2 += unbiased(n) ? 1 : 0;
}
System.out.format("%d: %2.2f%% %2.2f%%\n",
n, 100.0*c1/M, 100.0*c2/M);
}
}
}</lang> Output:
3: 33,11% 50,23% 4: 24.97% 49.78% 5: 20.05% 50.00% 6: 17.00% 49.88%
Liberty BASIC
<lang lb> for N =3 to 6 ' bias as defined
tests =1E5 ' number of tests to do
print " Biased bit-string, '1' chosen on average once out of "; N; " times . . . "
countZeros =0: countOnes =0
for j =1 to tests
b =randN( N)
if b =1 then countOnes =countOnes +1 else countZeros =countZeros +1
next j
print " "; countZeros; " zeros & "; countOnes; " ones. Ratio ="; countOnes /tests
print " Unbiased bit-string . . . "
countZeros =0: countOnes =0
for j =1 to tests
b =unBiased( N)
if b =1 then countOnes =countOnes +1 else countZeros =countZeros +1
next j
print " "; countZeros; " zeros & "; countOnes; " ones. Ratio ="; countOnes /tests print
next N
print " DONE."
end ' _____________________________________________________
function randN( n)
if rnd( 1) <( 1 /n) then randN =1 else randN =0
end function
function unBiased( n)
do
n1 =randN( n)
n2 =randN( n)
loop until n1 <>n2
unBiased =n1
end function </lang>
Output:
Biased bit-string, '1' chosen once out of 3 times . . . 664236 zeros & 335764 ones. Ratio =0.335764 Unbiased bit-string . . . 500349 zeros & 499651 ones. Ratio =0.499651 Biased bit-string, '1' chosen once out of 4 times . . . 748122 zeros & 251878 ones. Ratio =0.251878 Unbiased bit-string . . . 499728 zeros & 500272 ones. Ratio =0.500272 Biased bit-string, '1' chosen once out of 5 times . . . 798517 zeros & 201483 ones. Ratio =0.201483 Unbiased bit-string . . . 500044 zeros & 499956 ones. Ratio =0.499956 Biased bit-string, '1' chosen once out of 6 times . . . 832096 zeros & 167904 ones. Ratio =0.167904 Unbiased bit-string . . . 500407 zeros & 499593 ones. Ratio =0.499593
Lua
<lang lua> local function randN(n)
return function() if math.random() < 1/n then return 1 else return 0 end end
end
local function unbiased(n)
local biased = randN (n)
return function()
local a, b = biased(), biased()
while a==b do
a, b = biased(), biased()
end
return a
end
end
local function demonstrate (samples)
for n = 3, 6 do
biased = randN(n)
unbias = unbiased(n)
local bcounts = {[0]=0,[1]=0}
local ucounts = {[0]=0,[1]=0}
for i=1, samples do
local bnum = biased()
local unum = unbias()
bcounts[bnum] = bcounts[bnum]+1
ucounts[unum] = ucounts[unum]+1
end
print(string.format("N = %d",n),
"# 0", "# 1",
"% 0", "% 1")
print("biased", bcounts[0], bcounts[1],
bcounts[0] / samples * 100,
bcounts[1] / samples * 100)
print("unbias", ucounts[0], ucounts[1],
ucounts[0] / samples * 100,
ucounts[1] / samples * 100)
end
end
demonstrate(100000) </lang>
Output:
N = 3 # 0 # 1 % 0 % 1 biased 66832 33168 66.832 33.168 unbias 50207 49793 50.207 49.793 N = 4 # 0 # 1 % 0 % 1 biased 75098 24902 75.098 24.902 unbias 49872 50128 49.872 50.128 N = 5 # 0 # 1 % 0 % 1 biased 80142 19858 80.142 19.858 unbias 50049 49951 50.049 49.951 N = 6 # 0 # 1 % 0 % 1 biased 83407 16593 83.407 16.593 unbias 49820 50180 49.82 50.18
Mathematica
<lang Mathematica>rand[bias_, n_] := 1 - Unitize@RandomInteger[bias - 1, n]
unbiased[bias_, n_] :=
DeleteCases[rand[bias, {n, 2}], {a_, a_}]All, 1</lang>
count = 1000000;
TableForm[
Table[{n, Total[rand[n, count]]/count // N,
Total[#]/Length[#] &@unbiased[n, count] // N}, {n, 3, 6}],
TableHeadings -> {None, {n, "biased", "unbiased"}}]
n biased unbiased
3 0.33312 0.500074
4 0.24932 0.499883
5 0.1998 0.498421
6 0.16620 0.49805
NetRexx
<lang NetRexx>/* NetRexx */ options replace format comments java crossref symbols binary
runSample(arg) return
-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ method biased(n = int) public static returns boolean
return Math.random() < 1.0 / n
-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ method unbiased(n = int) public static returns boolean
a = boolean b = boolean loop until a \= b a = biased(n) b = biased(n) end return a
-- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ method runSample(arg) private static
parse arg Mx .
if Mx.length <= 0 then Mx = 50000
M = int Mx
loop n = int 3 to 6
c1 = int 0
c2 = int 0
loop for M
if biased(n) then c1 = c1 + 1
if unbiased(n) then c2 = c2 + 1
end
say Rexx(n).right(3)':' Rexx(100.0 * c1 / M).format(6, 2)'%' Rexx(100.0 * c2 / M).format(6, 2)'%'
end n
return
</lang> Output:
3: 32.78% 49.98% 4: 24.72% 50.31% 5: 19.95% 50.34% 6: 17.20% 50.20%
Nim
<lang nim>import math, strutils randomize()
template newSeqWith(len: int, init: expr): expr =
var result {.gensym.} = newSeq[type(init)](len)
for i in 0 .. <len:
result[i] = init
result
proc randN(n): (proc: range[0..1]) =
result = proc(): range[0..1] = if random(n) == 0: 1 else: 0
proc unbiased(biased): range[0..1] =
var (this, that) = (biased(), biased()) while this == that: this = biased() that = biased() return this
for n in 3..6:
var biased = randN(n)
var v = newSeqWith(1_000_000, biased())
var cnt0, cnt1 = 0
for x in v:
if x == 0: inc cnt0
else: inc cnt1
echo "Biased(",n,") = count1=",cnt1,", count0=",cnt0,", percent=",
formatFloat(100 * float(cnt1)/float(cnt1+cnt0), ffDecimal, 3)
v = newSeqWith(1_000_000, unbiased(biased))
cnt0 = 0
cnt1 = 0
for x in v:
if x == 0: inc cnt0
else: inc cnt1
echo " Unbiased = count1=",cnt1,", count0=",cnt0,", percent=",
formatFloat(100 * float(cnt1)/float(cnt1+cnt0), ffDecimal, 3)</lang>
Output:
Biased(3) = count1=332805, count0=667195, percent=33.281 Unbiased = count1=500157, count0=499843, percent=50.016 Biased(4) = count1=249575, count0=750425, percent=24.957 Unbiased = count1=500072, count0=499928, percent=50.007 Biased(5) = count1=199537, count0=800463, percent=19.954 Unbiased = count1=499396, count0=500604, percent=49.940 Biased(6) = count1=166728, count0=833272, percent=16.673 Unbiased = count1=499712, count0=500288, percent=49.971
OCaml
<lang ocaml>let randN n =
if Random.int n = 0 then 1 else 0
let rec unbiased n =
let a = randN n in if a <> randN n then a else unbiased n
let () =
Random.self_init();
let n = 50_000 in
for b = 3 to 6 do
let cb = ref 0 in
let cu = ref 0 in
for i = 1 to n do
cb := !cb + (randN b);
cu := !cu + (unbiased b);
done;
Printf.printf "%d: %5.2f%% %5.2f%%\n"
b (100.0 *. float !cb /. float n) (100.0 *. float !cu /. float n)
done</lang>
Output:
3: 33.07% 49.90% 4: 25.11% 49.85% 5: 19.82% 50.09% 6: 16.51% 50.51%
PARI/GP
GP's random number generation is high-quality, using Brent's XORGEN. Thus this program is slow: the required 400,000 unbiased numbers generated through this bias/unbias scheme take nearly a second. This requires about two million calls to random, which in turn generate a total of about three million calls to the underlying random number generator through the rejection strategy. The overall efficiency of the scheme is 0.8% for 32-bit and 0.4% for 64-bit...
<lang parigp>randN(N)=!random(N);
unbiased(N)={
my(a,b); while(1, a=randN(N); b=randN(N); if(a!=b, return(a)) )
}; for(n=3,6,print(n"\t"sum(k=1,1e5,unbiased(n))"\t"sum(k=1,1e5,randN(n))))</lang>
Output:
3 49997 33540 4 49988 24714 5 50143 20057 6 49913 16770
Perl
<lang perl>sub randn {
my $n = shift;
return int(rand($n) / ($n - 1));
}
for my $n (3 .. 6) {
print "Bias $n: ";
my (@raw, @fixed);
for (1 .. 10000) {
my $x = randn($n);
$raw[$x]++;
$fixed[$x]++ if randn($n) != $x
}
print "@raw, ";
printf("%3g+-%.3g%%\tfixed: ", $raw[0]/100,
100 * sqrt($raw[0] * $raw[1]) / ($raw[0] + $raw[1])**1.5);
print "@fixed, ";
printf("%3g+-%.3g%%\n", 100*$fixed[0]/($fixed[0] + $fixed[1]),
100 * sqrt($fixed[0] * $fixed[1]) / ($fixed[0] + $fixed[1])**1.5);
}</lang> Output:
Bias 3: 6684 3316, 66.84+-0.471% fixed: 2188 2228, 49.5471+-0.752% Bias 4: 7537 2463, 75.37+-0.431% fixed: 1924 1845, 51.048+-0.814% Bias 5: 7993 2007, 79.93+-0.401% fixed: 1564 1597, 49.478+-0.889% Bias 6: 8309 1691, 83.09+-0.375% fixed: 1403 1410, 49.8756+-0.943%
Perl 6
<lang perl6>sub randN ( $n where 3..6 ) {
return ( $n.rand / ($n - 1) ).Int;
}
sub unbiased ( $n where 3..6 ) {
my $n1;
repeat { $n1 = randN($n) } until $n1 != randN($n);
return $n1;
}
my $iterations = 1000; for 3 .. 6 -> $n {
my ( @raw, @fixed );
for ^$iterations {
@raw[ randN($n) ]++;
@fixed[ unbiased($n) ]++;
}
printf "N=%d randN: %s, %4.1f%% unbiased: %s, %4.1f%%\n",
$n, map { .perl, .[1] * 100 / $iterations }, @raw, @fixed;
}</lang>
Output:
N=3 randN: [676, 324], 32.4% unbiased: [517, 483], 48.3% N=4 randN: [734, 266], 26.6% unbiased: [489, 511], 51.1% N=5 randN: [792, 208], 20.8% unbiased: [494, 506], 50.6% N=6 randN: [834, 166], 16.6% unbiased: [514, 486], 48.6%
Phix
Copy of Euphoria <lang Phix>function randN(integer N)
return rand(N) = 1
end function
function unbiased(integer N) integer a
while 1 do
a = randN(N)
if a!=randN(N) then
return a
end if
end while
end function
constant n = 10000 integer cb, cu for b=3 to 6 do
cb = 0
cu = 0
for i=1 to n do
cb += randN(b)
cu += unbiased(b)
end for
printf(1, "%d: %5.2f%% %5.2f%%\n", {b, 100 * cb / n, 100 * cu / n})
end for</lang>
- Output:
3: 32.83% 50.34% 4: 24.78% 50.01% 5: 20.21% 49.71% 6: 16.68% 49.67%
PicoLisp
<lang PicoLisp>(de randN (N)
(if (= 1 (rand 1 N)) 1 0) )
(de unbiased (N)
(use (A B)
(while
(=
(setq A (randN N))
(setq B (randN N)) ) )
A ) )</lang>
Test: <lang PicoLisp>(for N (range 3 6)
(tab (2 1 7 2 7 2)
N ":"
(format
(let S 0 (do 10000 (inc 'S (randN N))))
2 )
"%"
(format
(let S 0 (do 10000 (inc 'S (unbiased N))))
2 )
"%" ) )</lang>
Output:
3: 33.21 % 50.48 % 4: 25.06 % 49.79 % 5: 20.04 % 49.75 % 6: 16.32 % 49.02 %
PL/I
<lang PL/I> test: procedure options (main); /* 20 Nov. 2012 */
randN: procedure(N) returns (bit (1));
declare N fixed (1);
declare random builtin;
declare r fixed (2) external initial (-1);
if r >= 0 then do; r = r-1; return ('0'b); end;
r = random()*2*N;
return ('1'b);
end randN;
random: procedure returns (bit(1));
declare (r1, r2) bit (1);
do until (r1 ^= r2);
r1 = randN(N); r2 = randN(N);
end;
return (r1);
end random;
declare (biasedrn, unbiasedrn) (100) bit (1); declare N fixed (1);
put ('N Biased Unbiased (tally of 100 random numbers)');
do N = 3 to 6;
biasedrn = randN(N); unbiasedrn = random;
put skip edit (N, sum(biasedrn), sum(unbiasedrn)) (F(1), 2 F(10));
end;
end test; </lang> Results:
N Biased Unbiased (tally of 100 random numbers) 3 24 42 4 18 47 5 16 41 6 11 53
PowerShell
<lang PowerShell> function randN ( [int]$N )
{
[int]( ( Get-Random -Maximum $N ) -eq 0 )
}
function unbiased ( [int]$N )
{
do {
$X = randN $N
$Y = randN $N
}
While ( $X -eq $Y )
return $X
}
</lang> Note: The [pscustomobject] type accelerator, used to simplify making the test output look pretty, requires version 3.0 or higher. <lang PowerShell> $Tests = 1000 ForEach ( $N in 3..6 )
{
$Biased = 0
$Unbiased = 0
ForEach ( $Test in 1..$Tests )
{
$Biased += randN $N
$Unbiased += unbiased $N
}
[pscustomobject]@{ N = $N
"Biased Ones out of $Test" = $Biased
"Unbiased Ones out of $Test" = $Unbiased }
}
</lang>
- Output:
N Biased Ones out of 1000 Unbiased Ones out of 1000 - ----------------------- ------------------------- 3 322 503 4 273 518 5 217 515 6 173 486
PureBasic
<lang PureBasic>Procedure biased(n)
If Random(n) <> 1 ProcedureReturn 0 EndIf ProcedureReturn 1
EndProcedure
Procedure unbiased(n)
Protected a, b Repeat a = biased(n) b = biased(n) Until a <> b ProcedureReturn a
EndProcedure
- count = 100000
Define n, m, output.s For n = 3 To 6
Dim b_count(1) Dim u_count(1) For m = 1 To #count x = biased(n) b_count(x) + 1 x = unbiased(n) u_count(x) + 1 Next output + "N = " + Str(n) + #LF$ output + " biased =>" + #tab$ + "#0=" + Str(b_count(0)) + #tab$ + "#1=" +Str(b_count(1)) output + #tab$ + " ratio=" + StrF(b_count(1) / #count * 100, 2) + "%" + #LF$ output + " unbiased =>" + #tab$ + "#0=" + Str(u_count(0)) + #tab$ + "#1=" + Str(u_count(1)) output + #tab$ + " ratio=" + StrF(u_count(1) / #count * 100, 2) + "%" + #LF$
Next MessageRequester("Biased and Unbiased random number results", output)</lang> Sample output:
--------------------------- Biased and Unbiased random number results --------------------------- N = 3 biased => #0=74856 #1=25144 ratio=25.14% unbiased => #0=50066 #1=49934 ratio=49.93% N = 4 biased => #0=80003 #1=19997 ratio=20.00% unbiased => #0=49819 #1=50181 ratio=50.18% N = 5 biased => #0=83256 #1=16744 ratio=16.74% unbiased => #0=50268 #1=49732 ratio=49.73% N = 6 biased => #0=85853 #1=14147 ratio=14.15% unbiased => #0=49967 #1=50033 ratio=50.03%
Python
<lang python>from __future__ import print_function import random
def randN(N):
" 1,0 random generator factory with 1 appearing 1/N'th of the time" return lambda: random.randrange(N) == 0
def unbiased(biased):
'uses a biased() generator of 1 or 0, to create an unbiased one'
this, that = biased(), biased()
while this == that: # Loop until 10 or 01
this, that = biased(), biased()
return this # return the first
if __name__ == '__main__':
from collections import namedtuple
Stats = namedtuple('Stats', 'count1 count0 percent')
for N in range(3, 7):
biased = randN(N)
v = [biased() for x in range(1000000)]
v1, v0 = v.count(1), v.count(0)
print ( "Biased(%i) = %r" % (N, Stats(v1, v0, 100. * v1/(v1 + v0))) )
v = [unbiased(biased) for x in range(1000000)]
v1, v0 = v.count(1), v.count(0)
print ( " Unbiased = %r" % (Stats(v1, v0, 100. * v1/(v1 + v0)), ) )</lang>
Sample output
Biased(3) = Stats(count1=331800, count0=668200, percent=33.18) Unbiased = Stats(count1=499740, count0=500260, percent=49.973999999999997) Biased(4) = Stats(count1=249770, count0=750230, percent=24.977) Unbiased = Stats(count1=499707, count0=500293, percent=49.970700000000001) Biased(5) = Stats(count1=199764, count0=800236, percent=19.976400000000002) Unbiased = Stats(count1=499456, count0=500544, percent=49.945599999999999) Biased(6) = Stats(count1=167561, count0=832439, percent=16.7561) Unbiased = Stats(count1=499963, count0=500037, percent=49.996299999999998)
R
<lang rsplus>randN = function(N) sample.int(N, 1) == 1
unbiased = function(f)
{while ((x <- f()) == f()) {}
x}
samples = 10000 print(t(round(d = 2, sapply(3:6, function(N) c(
N = N, biased = mean(replicate(samples, randN(N))), unbiased = mean(replicate(samples, unbiased(function() randN(N)))))))))</lang>
Sample output:
N biased unbiased [1,] 3 0.32 0.50 [2,] 4 0.24 0.50 [3,] 5 0.21 0.49 [4,] 6 0.16 0.51
Racket
<lang racket>
- lang racket
- Using boolean #t/#f instead of 1/0
(define ((randN n)) (zero? (random n))) (define ((unbiased biased))
(let loop () (let ([r (biased)]) (if (eq? r (biased)) (loop) r))))
- Counts
(define N 1000000) (for ([n (in-range 3 7)])
(define (try% R) (round (/ (for/sum ([i N]) (if (R) 1 0)) N 1/100)))
(define biased (randN n))
(printf "Count: ~a => Biased: ~a%; Unbiased: ~a%.\n"
n (try% biased) (try% (unbiased biased))))
</lang>
- Output:
Count: 3 => Biased: 33%; Unbiased: 50%. Count: 4 => Biased: 25%; Unbiased: 50%. Count: 5 => Biased: 20%; Unbiased: 50%. Count: 6 => Biased: 17%; Unbiased: 50%.
REXX
<lang rexx>/*REXX program generates unbiased random numbers and displays the results to terminal.*/ parse arg # R seed . /*get optional parameters from the CL. */ if #== | #=="," then #=1000 /*# the number of SAMPLES to be used.*/ if R== | R=="," then R=6 /*R the high number for the range. */ if datatype(seed, 'W') then call random ,,seed /*Not specified? Use for RANDOM seed. */ w=12; pad=left(,5) /*width of columnar output; indentation*/ dash='─'; @b="biased"; @ub='un'@b /*literals for the SAY column headers. */ say pad c('N',5) c(@b) c(@b'%') c(@ub) c(@ub"%") c('samples') /*six column header.*/ dash=
do N=3 to R; b=0; u=0; do j=1 for #; b=b+randN(N)
u=u+unbiased()
end /*j*/
say pad c(N,5) c(b) pct(b) c(u) pct(u) c(#)
end /*N*/
exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ c: return center( arg(1), word(arg(2) w, 1), left(dash, 1) ) pct: return c( format(arg(1) / # * 100, , 2)'%' ) /*two decimal digits.*/ randN: parse arg z; return random(1, z)==z unbiased: do until x\==randN(N); x=randN(N); end /*until*/; return x</lang> output when using the default inputs:
──N── ───biased─── ──biased%─── ──unbiased── ─unbiased%── ──samples───
3 348 34.80% 541 54.10% 1000
4 259 25.90% 479 47.90% 1000
5 188 18.80% 475 47.50% 1000
6 178 17.80% 488 48.80% 1000
output when using the input of: 10000
──N── ───biased─── ──biased%─── ──unbiased── ─unbiased%── ──samples───
3 3435 34.35% 4995 49.95% 10000
4 2535 25.35% 4957 49.57% 10000
5 2019 20.19% 4958 49.58% 10000
6 1644 16.44% 4982 49.82% 10000
output when using the input of: 100000 30
──N── ───biased─── ──biased%─── ──unbiased── ─unbiased%── ──samples───
3 33301 33.30% 50066 50.07% 100000
4 25359 25.36% 49401 49.40% 100000
5 20026 20.03% 49966 49.97% 100000
6 16579 16.58% 49956 49.96% 100000
7 14294 14.29% 50008 50.01% 100000
8 12402 12.40% 50479 50.48% 100000
9 11138 11.14% 50099 50.10% 100000
10 9973 9.97% 49988 49.99% 100000
11 9062 9.06% 50009 50.01% 100000
12 8270 8.27% 49929 49.93% 100000
13 7704 7.70% 49876 49.88% 100000
14 7223 7.22% 50414 50.41% 100000
15 6725 6.73% 50043 50.04% 100000
16 6348 6.35% 50252 50.25% 100000
17 5900 5.90% 49977 49.98% 100000
18 5583 5.58% 49991 49.99% 100000
19 5139 5.14% 49958 49.96% 100000
20 4913 4.91% 50198 50.20% 100000
21 4714 4.71% 49892 49.89% 100000
22 4517 4.52% 49760 49.76% 100000
23 4226 4.23% 50021 50.02% 100000
24 4174 4.17% 50141 50.14% 100000
25 4005 4.01% 49816 49.82% 100000
26 3890 3.89% 49819 49.82% 100000
27 3705 3.71% 50036 50.04% 100000
28 3567 3.57% 49665 49.67% 100000
29 3481 3.48% 50094 50.09% 100000
30 3355 3.36% 49831 49.83% 100000
Ring
<lang ring> for n = 3 to 6
biased = 0
unb = 0
for i = 1 to 10000
biased += randN(n)
unb += unbiased(n)
next
see "N = " + n + " : biased = " + biased/100 + "%, unbiased = " + unb/100 + "%" + nl
next
func unbiased nr
while 1
a = randN(nr)
if a != randN(nr) return a ok
end
func randN m
m = (random(m) = 1)
return m
</lang> Output:
N = 3 : biased = 25.38%, unbiased = 50.12% N = 4 : biased = 20.34%, unbiased = 49.17% N = 5 : biased = 16.65%, unbiased = 48.86% N = 6 : biased = 13.31%, unbiased = 49.96%
Ruby
<lang ruby>def rand_n(bias)
rand(bias) == 0 ? 1 : 0
end
def unbiased(bias)
a, b = rand_n(bias), rand_n(bias) until a != b #loop until a and b are 0,1 or 1,0 a
end
runs = 1_000_000 keys = %i(bias biased unbiased) #use [:bias,:biased,:unbiased] in Ruby < 2.0 puts keys.join("\t")
(3..6).each do |bias|
counter = Hash.new(0) # counter will respond with 0 when key is not known
runs.times do
counter[:biased] += 1 if rand_n(bias) == 1 #the first time, counter has no key for :biased, so it will respond 0
counter[:unbiased] += 1 if unbiased(bias) == 1
end
counter[:bias] = bias
puts counter.values_at(*keys).join("\t")
end</lang>
- Output:
bias biased unbiased 3 333043 500161 4 249133 499393 5 199767 500354 6 166163 499809
Scala
<lang scala>def biased( n:Int ) = scala.util.Random.nextFloat < 1.0 / n
def unbiased( n:Int ) = { def loop : Boolean = { val a = biased(n); if( a != biased(n) ) a else loop }; loop }
for( i <- (3 until 7) ) println {
val m = 50000
var c1,c2 = 0
(0 until m) foreach { j => if( biased(i) ) c1 += 1; if( unbiased(i) ) c2 += 1 }
"%d: %2.2f%% %2.2f%%".format(i, 100.0*c1/m, 100.0*c2/m)
}</lang>
- Output:
3: 33.09% 49.79% 4: 24.92% 49.92% 5: 19.75% 49.92% 6: 16.67% 50.23%
Seed7
<lang seed7>$ include "seed7_05.s7i";
include "float.s7i";
const func integer: randN (in integer: n) is
return ord(rand(1, n) = 1);
const func integer: unbiased (in integer: n) is func
result
var integer: unbiased is 0;
begin
repeat
unbiased := randN(n);
until unbiased <> randN(n);
end func;
const proc: main is func
local
const integer: tests is 50000;
var integer: n is 0;
var integer: sumBiased is 0;
var integer: sumUnbiased is 0;
var integer: count is 0;
begin
for n range 3 to 6 do
sumBiased := 0;
sumUnbiased := 0;
for count range 1 to tests do
sumBiased +:= randN(n);
sumUnbiased +:= unbiased(n);
end for;
writeln(n <& ": " <& flt(100 * sumBiased) / flt(tests) digits 3 lpad 6 <&
" " <& flt(100 * sumUnbiased) / flt(tests) digits 3 lpad 6);
end for;
end func;</lang>
Output:
3: 33.004 50.024 4: 25.158 50.278 5: 20.186 49.978 6: 16.570 49.936
Sidef
<lang ruby>func randN (n) {
n.rand / (n-1) -> int
}
func unbiased(n) {
var n1 = nil
do { n1 = randN(n) } while (n1 == randN(n))
return n1
}
var iterations = 1000
for n in (3..6) {
var raw = []
var fixed = []
iterations.times {
raw[ randN(n) ] := 0 ++
fixed[ unbiased(n) ] := 0 ++
}
printf("N=%d randN: %s, %4.1f%% unbiased: %s, %4.1f%%\n",
n, [raw, fixed].map {|a| (a.dump, a[1] * 100 / iterations) }...)
}</lang>
- Output:
N=3 randN: [661, 339], 33.9% unbiased: [497, 503], 50.3% N=4 randN: [765, 235], 23.5% unbiased: [493, 507], 50.7% N=5 randN: [812, 188], 18.8% unbiased: [509, 491], 49.1% N=6 randN: [820, 180], 18.0% unbiased: [510, 490], 49.0%
Tcl
<lang tcl># 1,0 random generator factory with 1 appearing 1/N'th of the time proc randN n {expr {rand()*$n < 1}}
- uses a biased generator of 1 or 0, to create an unbiased one
proc unbiased {biased} {
while 1 {
if {[set a [eval $biased]] != [eval $biased]} {return $a}
}
}
for {set n 3} {$n <= 6} {incr n} {
set biased [list randN $n]
for {set i 0;array set c {0 0 1 0}} {$i < 1000000} {incr i} {
incr c([eval $biased])
} puts [format " biased %d => #0=%d #1=%d ratio=%.2f%%" $n $c(0) $c(1) \
[expr {100.*$c(1)/$i}]]
for {set i 0;array set c {0 0 1 0}} {$i < 1000000} {incr i} {
incr c([unbiased $biased])
} puts [format "unbiased %d => #0=%d #1=%d ratio=%.2f%%" $n $c(0) $c(1) \
[expr {100.*$c(1)/$i}]] }</lang> Sample output:
biased 3 => #0=667076 #1=332924 ratio=33.29% unbiased 3 => #0=500263 #1=499737 ratio=49.97% biased 4 => #0=750470 #1=249530 ratio=24.95% unbiased 4 => #0=500644 #1=499356 ratio=49.94% biased 5 => #0=800243 #1=199757 ratio=19.98% unbiased 5 => #0=500878 #1=499122 ratio=49.91% biased 6 => #0=833623 #1=166377 ratio=16.64% unbiased 6 => #0=500518 #1=499482 ratio=49.95%
zkl
<lang zkl>fcn randN(N){ (not (0).random(N)).toInt() } fcn unbiased(randN){ while((a:=randN())==randN()){} a }</lang> <lang zkl>const Z=0d100_000; foreach N in ([3..6]){
"%d: biased: %3.2f%%, unbiased: %3.2f%%".fmt(N,
(0).reduce(Z,'wrap(s,_){ s+randN(N) },0.0)/Z*100,
(0).reduce(Z,'wrap(s,_){ s+unbiased(randN.fp(N)) },0.0)/Z*100)
.println();
}</lang>
- Output:
3: biased: 33.46%, unbiased: 49.80% 4: biased: 24.95%, unbiased: 50.01% 5: biased: 19.89%, unbiased: 50.18% 6: biased: 16.75%, unbiased: 50.22%
- Programming Tasks
- Solutions by Programming Task
- Ada
- Aime
- AutoHotkey
- Examples needing attention
- BBC BASIC
- C
- C sharp
- Clojure
- CoffeeScript
- Common Lisp
- D
- Elixir
- ERRE
- Euphoria
- F Sharp
- Fortran
- GAP
- Go
- Haskell
- Icon
- Unicon
- J
- Java
- Liberty BASIC
- Lua
- Mathematica
- NetRexx
- Nim
- OCaml
- PARI/GP
- Perl
- Perl 6
- Phix
- PicoLisp
- PL/I
- PowerShell
- PureBasic
- Python
- R
- Racket
- REXX
- Ring
- Ruby
- Scala
- Seed7
- Sidef
- Tcl
- Zkl
- GUISS/Omit
- Randomness