Unbias a random generator
(Redirected from Unbalance a random generator)
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 occurring, , is not the same as , the probability of a zero occurring, the probability of the occurrence of a one followed by a zero is × , assuming independence. 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.
11l
F randN(n)
‘1,0 random generator factory with 1 appearing 1/n'th of the time’
R () -> random:(@=n) == 0
F unbiased(biased)
‘uses a biased() generator of 1 or 0, to create an unbiased one’
V (this, that) = (biased(), biased())
L this == that
(this, that) = (biased(), biased())
R this
L(n) 3..6
V biased = randN(n)
V v = (0.<1000000).map(x -> @biased())
V (v1, v0) = (v.count(1), v.count(0))
print(‘Biased(#.): count1=#., count0=#., percent=#.2’.format(n, v1, v0, 100.0 * v1 / (v1 + v0)))
v = (0.<1000000).map(x -> unbiased(@biased))
(v1, v0) = (v.count(1), v.count(0))
print(‘ Unbiased: count1=#., count0=#., percent=#.2’.format(v1, v0, 100.0 * v1 / (v1 + v0)))- Output:
Biased(3): count1=332946, count0=667054, percent=33.29 Unbiased: count1=499751, count0=500249, percent=49.98 Biased(4): count1=250561, count0=749439, percent=25.06 Unbiased: count1=500576, count0=499424, percent=50.06 Biased(5): count1=200056, count0=799944, percent=20.01 Unbiased: count1=499975, count0=500025, percent=50.00 Biased(6): count1=165953, count0=834047, percent=16.60 Unbiased: count1=500104, count0=499896, percent=50.01
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;
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
integer
biased(integer bias)
{
1 ^ min(drand(bias - 1), 1);
}
integer
unbiased(integer bias)
{
integer a;
while ((a = biased(bias)) == biased(bias)) {
}
a;
}
integer
main(void)
{
integer b, n, cb, cu, i;
n = 10000;
b = 3;
while (b <= 6) {
i = cb = cu = 0;
while ((i += 1) <= n) {
cb += biased(b);
cu += unbiased(b);
}
o_form("bias ~: /d2p2/%% vs /d2p2/%%\n", b, 100r * cb / n,
100r * cu / n);
b += 1;
}
0;
}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
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
BASIC
BASIC256
function randN (n)
if int(rand * n) + 1 <> 1 then return 0 else return 1
end function
function unbiased (n)
do
a = randN (n)
b = randN (n)
until a <> b
return a
end function
numveces = 100000
print "Resultados de números aleatorios sesgados e imparciales" + chr(10)
for n = 3 to 6
dim b_numveces(n) fill 0
dim u_numveces(n) fill 0
for m = 1 to numveces
x = randN (n)
b_numveces[x] += 1
x = unbiased (n)
u_numveces[x] += 1
next m
print "N = "; n
print " Biased =>", "#0="; (b_numveces[0]); " #1="; (b_numveces[1]); " ratio = "; (b_numveces[1]/numveces*100); "%"
print "Unbiased =>", "#0="; (u_numveces[0]); " #1="; (u_numveces[1]); " ratio = "; (u_numveces[1]/numveces*100); "%"
next n
end- Output:
Resultados de números aleatorios sesgados e imparciales N = 3 Biased => #0=66625 #1=33375 ratio = 33.375% Unbiased => #0=50026 #1=49974 ratio = 49.974% N = 4 Biased => #0=74988 #1=25012 ratio = 25.012% Unbiased => #0=49809 #1=50191 ratio = 50.191% N = 5 Biased => #0=79893 #1=20107 ratio = 20.107% Unbiased => #0=50102 #1=49898 ratio = 49.898% N = 6 Biased => #0=83432 #1=16568 ratio = 16.568% Unbiased => #0=50091 #1=49909 ratio = 49.909%
BBC BASIC
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)
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
#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;
}
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++
#include <iostream>
#include <random>
std::default_random_engine generator;
bool biased(int n) {
std::uniform_int_distribution<int> distribution(1, n);
return distribution(generator) == 1;
}
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 = biased(n);
flip2 = biased(n);
} while (flip1 == flip2);
return flip1;
}
int main() {
for (size_t n = 3; n <= 6; n++) {
int biasedZero = 0;
int biasedOne = 0;
int unbiasedZero = 0;
int unbiasedOne = 0;
for (size_t i = 0; i < 100000; i++) {
if (biased(n)) {
biasedOne++;
} else {
biasedZero++;
}
if (unbiased(n)) {
unbiasedOne++;
} else {
unbiasedZero++;
}
}
std::cout << "(N = " << n << ")\n";
std::cout << "Biased:\n";
std::cout << " 0 = " << biasedZero << "; " << biasedZero / 1000.0 << "%\n";
std::cout << " 1 = " << biasedOne << "; " << biasedOne / 1000.0 << "%\n";
std::cout << "Unbiased:\n";
std::cout << " 0 = " << unbiasedZero << "; " << unbiasedZero / 1000.0 << "%\n";
std::cout << " 1 = " << unbiasedOne << "; " << unbiasedOne / 1000.0 << "%\n";
std::cout << '\n';
}
return 0;
}
- Output:
(N = 3) Biased: 0 = 66614; 66.614% 1 = 33386; 33.386% Unbiased: 0 = 49965; 49.965% 1 = 50035; 50.035% (N = 4) Biased: 0 = 75032; 75.032% 1 = 24968; 24.968% Unbiased: 0 = 50030; 50.03% 1 = 49970; 49.97% (N = 5) Biased: 0 = 80178; 80.178% 1 = 19822; 19.822% Unbiased: 0 = 49878; 49.878% 1 = 50122; 50.122% (N = 6) Biased: 0 = 83494; 83.494% 1 = 16506; 16.506% Unbiased: 0 = 50085; 50.085% 1 = 49915; 49.915%
C#
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;
}
}
}
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
(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])
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
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
(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))))
output
3: unbiased 4992 biased 3361 4: unbiased 4988 biased 2472 5: unbiased 5019 biased 1987 6: unbiased 4913 biased 1658
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);
}
- Output:
3: 33.441% 49.964% 4: 24.953% 49.910% 5: 19.958% 49.987% 6: 16.660% 49.890%
EasyLang
func biased n .
return if randomf < 1 / n
.
func unbiased n .
repeat
a = biased n
b = biased n
until a <> b
.
return a
.
m = 50000
for n = 3 to 6
c1 = 0
c2 = 0
for i to m
c1 += biased n
c2 += unbiased n
.
print n & ": " & 100 * c1 / m & " " & 100 * c2 / m
.
- Output:
3: 33.10 50.14 4: 25.03 50.14 5: 19.81 50.04 6: 16.75 49.83
Elena
ELENA 6.x :
import extensions;
extension op : IntNumber
{
bool randN()
= randomGenerator.nextInt(self) == 0;
get bool Unbiased()
{
bool flip1 := self.randN();
bool flip2 := self.randN();
while (flip1 == flip2)
{
flip1 := self.randN();
flip2 := self.randN()
};
^ flip1
}
}
public program()
{
for(int n := 3; n <= 6; n += 1)
{
int biasedZero := 0;
int biasedOne := 0;
int unbiasedZero := 0;
int unbiasedOne := 0;
for(int i := 0; i < 100000; i += 1)
{
if(n.randN()) { biasedOne += 1 } else { biasedZero += 1 };
if(n.Unbiased) { unbiasedOne += 1 } else { unbiasedZero += 1 }
};
console
.printLineFormatted("(N = {0}):".padRight(17) + "# of 0"$9"# of 1"$9"% of 0"$9"% of 1", n)
.printLineFormatted("Biased:".padRight(15) + "{0}"$9"{1}"$9"{2}"$9"{3}",
biasedZero, biasedOne, biasedZero / 1000, biasedOne / 1000)
.printLineFormatted("Unbiased:".padRight(15) + "{0}"$9"{1}"$9"{2}"$9"{3}",
unbiasedZero, unbiasedOne, unbiasedZero / 1000, unbiasedOne / 1000)
}
}- Output:
(N = 3): # of 0 # of 1 % of 0 % of 1 Biased: 66793 33207 66 33 Unbiased: 49965 50035 49 50 (N = 4): # of 0 # of 1 % of 0 % of 1 Biased: 75233 24767 75 24 Unbiased: 50106 49894 50 49 (N = 5): # of 0 # of 1 % of 0 % of 1 Biased: 80209 19791 80 19 Unbiased: 50080 49920 50 49 (N = 6): # of 0 # of 1 % of 0 % of 1 Biased: 83349 16651 83 16 Unbiased: 49699 50301 49 50
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
- 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
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- 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
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 forOutput:
3: 33.68% 49.94% 4: 24.93% 50.48% 5: 20.32% 49.97% 6: 16.98% 50.05%
F#
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
- Output:
3: 33.26% 49.97% 4: 25.02% 50.22% 5: 19.98% 50.00% 6: 16.64% 49.69%
Factor
USING: formatting kernel math math.ranges random sequences ;
IN: rosetta-code.unbias
: randN ( n -- m ) random zero? 1 0 ? ;
: unbiased ( n -- m )
dup [ randN ] dup bi 2dup = not
[ drop nip ] [ 2drop unbiased ] if ;
: test-generator ( quot -- x )
[ 1,000,000 dup ] dip replicate sum 100 * swap / ; inline
: main ( -- )
3 6 [a,b] [
dup [ randN ] [ unbiased ] bi-curry
[ test-generator ] bi@ "%d: %.2f%% %.2f%%\n" printf
] each ;
MAIN: main
- Output:
3: 33.25% 50.03% 4: 24.98% 50.02% 5: 20.03% 50.04% 6: 16.66% 49.99%
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
Output:
3: 33.337% 49.971% 4: 24.945% 49.944% 5: 19.971% 49.987% 6: 16.688% 50.097%
FreeBASIC
Function randN (n As Ubyte) As Ubyte
If Int(Rnd * n) + 1 <> 1 Then Return 0 Else Return 1
End Function
Function unbiased (n As Ubyte) As Ubyte
Dim As Ubyte a, b
Do
a = randN (n)
b = randN (n)
Loop Until a <> b
Return a
End Function
Const count = 100000
Dim x As Ubyte
Randomize Timer
Print "Resultados de n";Chr(163);!"meros aleatorios sesgados e imparciales\n"
For n As Ubyte = 3 To 6
Dim As Integer b_count(1)
Dim As Integer u_count(1)
For m As Integer = 1 To count
x = randN (n)
b_count(x) += 1
x = unbiased (n)
u_count(x) += 1
Next m
Print "N ="; n
Print " Biased =>", "#0="; Str(b_count(0)), "#1="; Str(b_count(1)),
Print Using "ratio = ##.##%"; (b_count(1) / count * 100)
Print "Unbiased =>", "#0="; Str(u_count(0)), "#1="; Str(u_count(1)),
Print Using "ratio = ##.##%"; (u_count(1) / count * 100)
Next n
Sleep
- Output:
Resultados de números aleatorios sesgados e imparciales N = 3 Biased => #0=66633 #1=33367 ratio = 33.37% Unbiased => #0=50029 #1=49971 ratio = 49.97% N = 4 Biased => #0=75068 #1=24932 ratio = 24.93% Unbiased => #0=50107 #1=49893 ratio = 49.89% N = 5 Biased => #0=79998 #1=20002 ratio = 20.00% Unbiased => #0=50049 #1=49951 ratio = 49.95% N = 6 Biased => #0=83195 #1=16805 ratio = 16.81% Unbiased => #0=50026 #1=49974 ratio = 49.97%
Fōrmulæ
Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text. Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation —i.e. XML, JSON— they are intended for storage and transfer purposes more than visualization and edition.
Programs in Fōrmulæ are created/edited online in its website.
In this page you can see and run the program(s) related to this task and their results. You can also change either the programs or the parameters they are called with, for experimentation, but remember that these programs were created with the main purpose of showing a clear solution of the task, and they generally lack any kind of validation.
Solution
Test cases
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 ] ]
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)
}
}
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:
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)]
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.
unbiased :: (MonadRandom m, Eq x) => m x -> m x
unbiased g = do x <- g
y <- g
if x /= y then return y else unbiased g
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 (unbiased (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:
main = forM_ [3..6] showCounts
where
showCounts b = do
r1 <- counts (randN b)
r2 <- counts (unbiased (randN b))
printf "n = %d biased: %d%% unbiased: %d%%\n" b r1 r2
counts g = (`div` 100) . length . filter (== 1) <$> replicateM 10000 g
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.
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
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
randN=: 0 = ?
unbiased=: i.@# { ::$: 2 | 0 3 -.~ _2 #.\ 4&* randN@# ]
Example use:
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
Some example counts (these are counts of the number of 1s which appear in a test involving 100 random numbers):
+/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
Note that these results are arbitrary. 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
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);
}
}
}
Output:
3: 33,11% 50,23% 4: 24.97% 49.78% 5: 20.05% 50.00% 6: 17.00% 49.88%
jq
Adapted from #Wren
Also works with gojq, the Go implementation of jq
In this entry, /dev/random is used as a source of entropy, and so a suitable invocation would be:
< /dev/random tr -cd '0-9' | fold -w 1 | jq -Mcnr -f unbias.jq
unbias.jq
### Utility Functions
# Output: a PRN in range(0;$n) where $n is .
def prn:
if . == 1 then 0
else . as $n
| (($n-1)|tostring|length) as $w
| [limit($w; inputs)] | join("") | tonumber
| if . < $n then . else ($n | prn) end
end;
def round: ((. * 100) | floor) / 100;
# input: n
# output: boolean, such that P(true) == 1/n
def biased:
prn == 1 | debug;
def unbiased:
. as $n
| {}
| until( .a != .b; {a: ($n|biased), b: ($n|biased)})
| .a;
def task($m):
def f(x;y;z): "\(x): \(y|round)% \(z|round)%";
range(3;7) as $n
| reduce range(0; $m) as $i ( {c1:0, c2:0};
if ($n|biased) then .c1 += 1 else . end
| if ($n|unbiased) then .c2 += 1 else . end)
| f($n; 100 * .c1 / $m; 100 * .c2 / $m);
task(50000)- Output:
3: 33.61% 50.09% 4: 25.48% 50.21% 5: 20.32% 50.47% 6: 16.35% 49.72%
Julia
using Printf
randN(N) = () -> rand(1:N) == 1 ? 1 : 0
function unbiased(biased::Function)
this, that = biased(), biased()
while this == that this, that = biased(), biased() end
return this
end
@printf "%2s | %10s | %5s | %5s | %8s" "N" "bias./unb." "1s" "0s" "pct ratio"
const nrep = 10000
for N in 3:6
biased = randN(N)
v = collect(biased() for __ in 1:nrep)
v1, v0 = count(v .== 1), count(v .== 0)
@printf("%2i | %10s | %5i | %5i | %5.2f%%\n", N, "biased", v1, v0, 100 * v1 / nrep)
v = collect(unbiased(biased) for __ in 1:nrep)
v1, v0 = count(v .== 1), count(v .== 0)
@printf("%2i | %10s | %5i | %5i | %5.2f%%\n", N, "unbiased", v1, v0, 100 * v1 / nrep)
end
- Output:
N | bias./unb. | 1s | 0s | pct ratio 3 | biased | 3286 | 6714 | 32.86% 3 | unbiased | 4986 | 5014 | 49.86% 4 | biased | 2473 | 7527 | 24.73% 4 | unbiased | 4986 | 5014 | 49.86% 5 | biased | 1992 | 8008 | 19.92% 5 | unbiased | 5121 | 4879 | 51.21% 6 | biased | 1663 | 8337 | 16.63% 6 | unbiased | 5040 | 4960 | 50.40%
Kotlin
// version 1.1.2
fun biased(n: Int) = Math.random() < 1.0 / n
fun unbiased(n: Int): Boolean {
var a: Boolean
var b: Boolean
do {
a = biased(n)
b = biased(n)
}
while (a == b)
return a
}
fun main(args: Array<String>) {
val m = 50_000
val f = "%d: %2.2f%% %2.2f%%"
for (n in 3..6) {
var c1 = 0
var c2 = 0
for (i in 0 until m) {
if (biased(n)) c1++
if (unbiased(n)) c2++
}
println(f.format(n, 100.0 * c1 / m, 100.0 * c2 / m))
}
}
Sample output:
3: 33.19% 50.19% 4: 25.29% 49.85% 5: 19.91% 50.07% 6: 16.71% 50.14%
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)
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 /Wolfram Language
rand[bias_, n_] := 1 - Unitize@RandomInteger[bias - 1, n]
unbiased[bias_, n_] := DeleteCases[rand[bias, {n, 2}], {a_, a_}][[All, 1]]
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"}}]
- Output:
n biased unbiased 3 0.33312 0.500074 4 0.24932 0.499883 5 0.1998 0.498421 6 0.16620 0.49805
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
Output:
3: 32.78% 49.98% 4: 24.72% 50.31% 5: 19.95% 50.34% 6: 17.20% 50.20%
Nim
import random, sequtils, strformat
type randProc = proc: range[0..1]
randomize()
proc randN(n: Positive): randProc =
result = proc: range[0..1] = ord(rand(n) == 0)
proc unbiased(biased: randProc): range[0..1] =
result = biased()
var that = biased()
while result == that:
result = biased()
that = biased()
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 = {100*cnt1 / (cnt1+cnt0):.3f}"
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 = {100*cnt1 / (cnt1+cnt0):.3f}"
- Output:
Biased(3) → count1 = 249654, count0 = 750346, percent = 24.965 Unbiased → count1 = 498631, count0 = 501369, percent = 49.863 Biased(4) → count1 = 200117, count0 = 799883, percent = 20.012 Unbiased → count1 = 499425, count0 = 500575, percent = 49.943 Biased(5) → count1 = 167218, count0 = 832782, percent = 16.722 Unbiased → count1 = 500020, count0 = 499980, percent = 50.002 Biased(6) → count1 = 143388, count0 = 856612, percent = 14.339 Unbiased → count1 = 500378, count0 = 499622, percent = 50.038
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
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...
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))))Output:
3 49997 33540 4 49988 24714 5 50143 20057 6 49913 16770
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);
}
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%
Phix
with javascript_semantics function randN(integer N) return rand(N) = 1 end function function unbiased(integer N) while true do integer a = randN(N) if a!=randN(N) then return a end if end while end function constant n = 100000 for b=3 to 6 do integer cb = 0, cu = 0 for i=1 to n do cb += randN(b) cu += unbiased(b) end for printf(1, "%d: biased:%.0f%% unbiased:%.0f%%\n", {b,(cb/n)*100,(cu/n)*100}) end for
- Output:
(Rounded to the nearest whole percent, of course)
3: biased:33% unbiased:50% 4: biased:25% unbiased:50% 5: biased:20% unbiased:50% 6: biased:17% unbiased:50%
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 ) )Test:
(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 )
"%" ) )Output:
3: 33.21 % 50.48 % 4: 25.06 % 49.79 % 5: 20.04 % 49.75 % 6: 16.32 % 49.02 %
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;Results:
N Biased Unbiased (tally of 100 random numbers) 3 24 42 4 18 47 5 16 41 6 11 53
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
}
Note: The [pscustomobject] type accelerator, used to simplify making the test output look pretty, requires version 3.0 or higher.
$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 }
}
- 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
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)
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
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)), ) )
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)
Quackery
$ "bigrat.qky" loadfile
[ random 0 = ] is randN ( n --> n )
[ dup randN
over randN
2dup = iff
2drop again
drop nip ] is unbias ( n --> n )
[ dup echo say " biased --> "
0
1000000 times
[ over randN if 1+ ]
nip 1000000 6 point$ echo$ ] is showbias ( n --> )
[ dup echo say " unbiased --> "
0
1000000 times
[ over unbias if 1+ ]
nip 1000000 6 point$ echo$ ] is showunbias ( n --> )
' [ 3 4 5 6 ]
witheach
[ dup cr
showbias cr
showunbias cr ]- Output:
3 biased --> 0.333225 3 unbiased --> 0.500147 4 biased --> 0.249658 4 unbiased --> 0.499851 5 biased --> 0.200169 5 unbiased --> 0.500073 6 biased --> 0.166804 6 unbiased --> 0.499045
R
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)))))))))
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
;; 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))))
- 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%.
Raku
(formerly Perl 6)
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 { .raku, .[1] * 100 / $iterations }, @raw, @fixed;
}
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%
REXX
/*REXX program generates unbiased random numbers and displays the results to terminal.*/
parse arg # R seed . /*obtain optional arguments 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 /*Specified? Then use for RANDOM seed.*/
dash='─'; @b="biased"; @ub='un'@b /*literals for the SAY column headers. */
say left('',5) ctr("N",5) ctr(@b) ctr(@b'%') ctr(@ub) ctr(@ub"%") ctr('samples')
dash=
do N=3 to R; b=0; u=0
do j=1 for #; b=b + randN(N); u=u + unbiased()
end /*j*/
say left('', 5) ctr(N, 5) ctr(b) pct(b) ctr(u) pct(u) ctr(#)
end /*N*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
ctr: return center( arg(1), word(arg(2) 12, 1), left(dash, 1)) /*show hdr│numbers.*/
pct: return ctr( format(arg(1) / # * 100, , 2)'%' ) /*2 decimal digits.*/
randN: parse arg z; return random(1, z)==z /*ret 1 if rand==Z.*/
unbiased: do until x\==randN(N); x=randN(N); end; return x /* " unbiased RAND*/
- 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
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 mOutput:
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%
RPL
≪ INV RAND ≥ ≫ ‘RandN’ STO
≪ 0 DUP WHILE DUP2 == REPEAT
DROP2 DUP RandN OVER RandN
END ROT DROP2
≫ ‘Unbiased’ STO
≪ 3 6 FOR n
(0,0) 1 10000 START
n RandN + n Unbiased i * →NUM +
NEXT 10000 / NEXT
≫ ‘SHOW’ STO
- Input:
SHOW
- Output:
4: (0.3272,0.4976) 3: (0.2504,0.4961) 2: (0.201,0.5008) 1: (0.166,0.4952)
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
- Output:
bias biased unbiased 3 333043 500161 4 249133 499393 5 199767 500354 6 166163 499809
Rust
#![feature(inclusive_range_syntax)]
extern crate rand;
use rand::Rng;
fn rand_n<R: Rng>(rng: &mut R, n: u32) -> usize {
rng.gen_weighted_bool(n) as usize // maps `false` to 0 and `true` to 1
}
fn unbiased<R: Rng>(rng: &mut R, n: u32) -> usize {
let mut bit = rand_n(rng, n);
while bit == rand_n(rng, n) {
bit = rand_n(rng, n);
}
bit
}
fn main() {
const SAMPLES: usize = 100_000;
let mut rng = rand::weak_rng();
println!(" Bias rand_n unbiased");
for n in 3..=6 {
let mut count_biased = 0;
let mut count_unbiased = 0;
for _ in 0..SAMPLES {
count_biased += rand_n(&mut rng, n);
count_unbiased += unbiased(&mut rng, n);
}
let b_percentage = 100.0 * count_biased as f64 / SAMPLES as f64;
let ub_percentage = 100.0 * count_unbiased as f64 / SAMPLES as f64;
println!(
"bias {}: {:0.2}% {:0.2}%",
n, b_percentage, ub_percentage
);
}
}
- Output:
Bias rand_n unbiased bias 3: 33.32% 49.80% bias 4: 25.22% 50.16% bias 5: 19.91% 50.00% bias 6: 16.66% 49.95%
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)
}
- Output:
3: 33.09% 49.79% 4: 24.92% 49.92% 5: 19.75% 49.92% 6: 16.67% 50.23%
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;Output:
3: 33.004 50.024 4: 25.158 50.278 5: 20.186 49.978 6: 16.570 49.936
Sidef
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) }...)
}
- 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
# 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}]]
}
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%
Visual Basic .NET
Module Module1
Dim random As New Random()
Function RandN(n As Integer) As Boolean
Return random.Next(0, n) = 0
End Function
Function Unbiased(n As Integer) As Boolean
Dim flip1 As Boolean
Dim flip2 As Boolean
Do
flip1 = RandN(n)
flip2 = RandN(n)
Loop While flip1 = flip2
Return flip1
End Function
Sub Main()
For n = 3 To 6
Dim biasedZero = 0
Dim biasedOne = 0
Dim unbiasedZero = 0
Dim unbiasedOne = 0
For i = 1 To 100000
If RandN(n) Then
biasedOne += 1
Else
biasedZero += 1
End If
If Unbiased(n) Then
unbiasedOne += 1
Else
unbiasedZero += 1
End If
Next
Console.WriteLine("(N = {0}):".PadRight(17) + "# of 0" + vbTab + "# of 1" + vbTab + "% of 0" + vbTab + "% of 1", n)
Console.WriteLine("(Biased: {0}):".PadRight(15) + "{0}" + vbTab + "{1}" + vbTab + "{2}" + vbTab + "{3}", biasedZero, biasedOne, biasedZero / 1000, biasedOne / 1000)
Console.WriteLine("(UnBiased: {0}):".PadRight(15) + "{0}" + vbTab + "{1}" + vbTab + "{2}" + vbTab + "{3}", unbiasedZero, unbiasedOne, unbiasedZero / 1000, unbiasedOne / 1000)
Next
End Sub
End Module
- Output:
(N = 3): # of 0 # of 1 % of 0 % of 1 (Biased: 66844): 66844 33156 66.844 33.156 (UnBiased: 50081):50081 49919 50.081 49.919 (N = 4): # of 0 # of 1 % of 0 % of 1 (Biased: 74827): 74827 25173 74.827 25.173 (UnBiased: 50095):50095 49905 50.095 49.905 (N = 5): # of 0 # of 1 % of 0 % of 1 (Biased: 79878): 79878 20122 79.878 20.122 (UnBiased: 50544):50544 49456 50.544 49.456 (N = 6): # of 0 # of 1 % of 0 % of 1 (Biased: 83449): 83449 16551 83.449 16.551 (UnBiased: 49907):49907 50093 49.907 50.093
Wren
import "random" for Random
import "./fmt" for Fmt
var rand = Random.new()
var biased = Fn.new { |n| rand.float() < 1 / n }
var unbiased = Fn.new { |n|
while (true) {
var a = biased.call(n)
var b = biased.call(n)
if (a != b) return a
}
}
var m = 50000
var f = "$d: $2.2f\% $2.2f\%"
for (n in 3..6) {
var c1 = 0
var c2 = 0
for (i in 0...m) {
if (biased.call(n)) c1 = c1 + 1
if (unbiased.call(n)) c2 = c2 + 1
}
Fmt.print(f, n, 100 * c1 / m, 100 * c2 / m)
}
- Output:
Sample output:
3: 33.36% 49.91% 4: 24.80% 49.89% 5: 20.07% 49.95% 6: 16.64% 50.17%
zkl
fcn randN(N){ (not (0).random(N)).toInt() }
fcn unbiased(randN){ while((a:=randN())==randN()){} a }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();
}- 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%
Categories:
- Programming Tasks
- Solutions by Programming Task
- 11l
- Ada
- Aime
- AutoHotkey
- Examples needing attention
- BASIC
- BASIC256
- BBC BASIC
- C
- C++
- C sharp
- Clojure
- CoffeeScript
- Common Lisp
- D
- EasyLang
- Elena
- Elixir
- ERRE
- Euphoria
- F Sharp
- Factor
- Fortran
- FreeBASIC
- Fōrmulæ
- GAP
- Go
- Haskell
- Icon
- Unicon
- J
- Java
- Jq
- Julia
- Kotlin
- Lua
- Mathematica
- Wolfram Language
- NetRexx
- Nim
- OCaml
- PARI/GP
- Perl
- Phix
- PicoLisp
- PL/I
- PowerShell
- PureBasic
- Python
- Quackery
- R
- Racket
- Raku
- REXX
- Ring
- RPL
- Ruby
- Rust
- Scala
- Seed7
- Sidef
- Tcl
- Visual Basic .NET
- Wren
- Wren-fmt
- Zkl
- GUISS/Omit
- Randomness