Unbias a random generator: Difference between revisions
m (output to pre) |
(→{{header|Mathematica}}: Marked incorrect. The task asks for specific output to be shown that is missing.) |
||
| Line 791: | Line 791: | ||
=={{header|Mathematica}}== |
=={{header|Mathematica}}== |
||
{{incorrect|Mathematica|The task asks for specific output to be shown that is missing.}} |
|||
<lang Mathematica>randN[n_] := If[Mod[#, 13] == 0, 1, 0] & /@ RandomInteger[13, n] |
<lang Mathematica>randN[n_] := If[Mod[#, 13] == 0, 1, 0] & /@ RandomInteger[13, n] |
||
Revision as of 22:33, 23 June 2012
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.
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
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>
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;
bool biased(in int n) /*nothrow*/ {
return uniform(0.0, 1.0) < (1.0 / n);
}
int unbiased(in int bias) /*nothrow*/ {
int a;
while ((a = biased(bias)) == biased(bias)) {}
return a;
}
void main() {
alias reduce!q{a + b} sum; /**/
enum int M = 500_000;
foreach (n; 3 .. 7) {
immutable a1 = sum(0, iota(M).map!(_=> biased(n))()),
a2 = sum(0, iota(M).map!(_=> unbiased(n))());
writefln("%d: %2.3f%% %2.3f%%", n,
100.0 * a1 / M, 100.0 * a2 / M);
}
}</lang>
- Output:
3: 33.441% 49.964% 4: 24.953% 49.910% 5: 19.958% 49.987% 6: 16.660% 49.890%
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%
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
Crappy implementation using IO
<lang haskell>import Control.Monad
import Random
import Data.IORef
import Text.Printf
randN :: Integer -> IO Bool randN n = randomRIO (1,n) >>= return . (== 1)
unbiased :: Integer -> IO Bool unbiased n = do
a <- randN n b <- randN n if a /= b then return a else unbiased n
main :: IO () main = forM_ [3..6] $ \n -> do
cb <- newIORef 0 cu <- newIORef 0 replicateM_ trials $ do b <- randN n u <- unbiased n when b $ modifyIORef cb (+ 1) when u $ modifyIORef cu (+ 1) tb <- readIORef cb tu <- readIORef cu printf "%d: %5.2f%% %5.2f%%\n" n (100 * fromIntegral tb / fromIntegral trials :: Double) (100 * fromIntegral tu / fromIntegral trials :: Double) where trials = 50000</lang>
Output:
3: 33.72% 50.08% 4: 25.26% 50.15% 5: 19.99% 50.07% 6: 16.67% 50.10%
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>randN[n_] := If[Mod[#, 13] == 0, 1, 0] & /@ RandomInteger[13, n]
unbiased[n_] := Module[{b = {}, a = 0},
While[Length[b] < n, a = randN[2];
If[Differences[a] != {0}, AppendTo[b, First[a]]]]; b
]</lang>
randN[10]
-> {0, 1, 0, 0, 0, 0, 0, 0, 0, 0}
unbiased[10]
->{0, 1, 0, 0, 0, 0, 1, 1, 1, 0}
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%
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 %
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)
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
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%
- Programming Tasks
- Solutions by Programming Task
- Ada
- AutoHotkey
- Examples needing attention
- C
- C sharp
- Clojure
- CoffeeScript
- Common Lisp
- D
- Euphoria
- Fortran
- GAP
- Go
- Haskell
- Icon
- Unicon
- J
- Java
- Liberty BASIC
- Lua
- Mathematica
- Mathematica examples needing attention
- OCaml
- PARI/GP
- Perl
- Perl 6
- PicoLisp
- PureBasic
- Python
- Seed7
- Tcl
- GUISS/Omit
- Randomness