Unbias a random generator: Difference between revisions

 

* 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 [http://en.wikipedia.org/wiki/Randomness_extractor#Von_Neumann_extractor Von Neumann debiasing], first described in a 1951 paper.

 

=={{header|Ada}}==

 

procedure Bias_Unbias is

Ada.Text_IO.New_Line;

end loop;

Output:<pre> I Biased% UnBiased%

3 32.87 49.80

=={{header|Aime}}==

{{trans|C}}

biased(integer bias)

{

}

 

}

 

}

 

b = 3;

while (b <= 6) {

cb += biased(b);

cu += unbiased(b);

}

 

o_form("bias ~: /d2p2/%% vs /d2p2/%%\n", b, 100r * cb / n,

 

b += 1;

}

 

Output:<pre>bias 3: 33.51% vs 50.27%

bias 4: 24.97% vs 49.99%

=={{header|AutoHotkey}}==

{{output?}}

Random, q, 0, 4

return q=4

MsgBox % "Unbiased probability of a 1 occurring: " Errorlevel/1000

StringReplace, junk, t, 1, , UseErrorLevel

 

=={{header|BBC BASIC}}==

biased% = 0

unbiased% = 0

= A%

Output:

<pre>

 

=={{header|C}}==

#include <stdlib.h>

 

return 0;

output

<pre>bias 3: 33.090% vs 49.710%

bias 5: 19.760% vs 49.650%

bias 6: 16.740% vs 50.030%</pre>

 

=={{header|C sharp}}==

 

 

namespace Unbias

}

}

 

'''Sample Output'''

 

=={{header|Clojure}}==

(if (< (rand 2) (/ n)) 0 1))

 

[4 0.87684 0.5023]

[5 0.90122 0.49728]

 

=={{header|CoffeeScript}}==

biased_rand_function = (n) ->

# return a function that returns 0/1 with

stats "biased", f_biased

stats "unbiased", f_unbiased

output

<pre>

 

=={{header|Common Lisp}}==

 

(defun unbiased (n)

(let ((u (loop repeat 10000 collect (unbiased n)))

(b (loop repeat 10000 collect (biased n))))

output

<pre>3: unbiased 4992 biased 3361

 

=={{header|D}}==

 

enum biased = (in int n) /*nothrow*/ => uniform01 < (1.0 / n);

M.iota.map!(_=> n.biased).sum * 100.0 / M,

M.iota.map!(_=> n.unbiased).sum * 100.0 / M);

{{out}}

<pre>3: 33.441% 49.964%

5: 19.958% 49.987%

6: 16.660% 49.890%</pre>

 

=={{header|Elixir}}==

def randN(n) do

end

def unbiased(n) do

end

end

IO.puts "N biased unbiased"

for n <- 3..6 do

 

{{out}}

5 0.2027 0.5041

6 0.1647 0.4912

</pre>

 

=={{header|Euphoria}}==

return rand(N) = 1

end function

end for

printf(1, "%d: %5.2f%% %5.2f%%\n", {b, 100 * cb / n, 100 * cu / n})

 

Output:

5: 20.32% 49.97%

6: 16.98% 50.05%

</pre>

 

=={{header|Fortran}}==

{{works with|Fortran|90 and later}}

implicit none

 

end function

 

Output:

<pre>3: 33.337% 49.971%

5: 19.971% 49.987%

6: 16.688% 50.097%</pre>

 

=={{header|GAP}}==

local v, rand;

v := [1 .. n - 1]*0;

# [ 4, 249851, 500663 ],

# [ 5, 200532, 500448 ],

=={{header|Go}}==

 

import (

u[1], u[0], float64(u[1])*100/samples)

}

Output:

<pre>

 

=={{header|Haskell}}==

import Text.Printf

 

 

 

 

Output:

 

<pre>

</pre>

 

This solution works in both languages. Both <tt>randN</tt> and

<tt>unbiased</tt> are generators in the Icon/Unicon sense.

iters := \A[1] | 10000

write("ratios of 0 to 1 from ",iters," trials:")

procedure randN(n)

repeat suspend if 1 = ?n then 1 else 0

and a sample run:

<pre>->ubrn 100000

 

=={{header|J}}==

 

Example use:

 

1 0 0 0 1 0 0 0 0 0

unbiased 10#6

 

Some example counts (these are counts of the number of 1s which appear in a test involving 100 random numbers):

 

30

+/randN 100#4

49

+/unbiased 100#6

 

 

=={{header|Java}}==

public static boolean biased(int n) {

return Math.random() < 1.0 / n;

}

}

Output:

<pre>3: 33,11% 50,23%

6: 17.00% 49.88%</pre>

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

</pre>

 

=={{header|Lua}}==

local function randN(n)

return function()

 

demonstrate(100000)

 

Output:

</pre>

 

TableForm[

Table[{n, Total[rand[n, count]]/count // N,

Total[#]/Length[#] &@unbiased[n, count] // N}, {n, 3, 6}],

3 0.33312 0.500074

4 0.24932 0.499883

5 0.1998 0.498421

 

=={{header|NetRexx}}==

{{trans|Java}}

options replace format comments java crossref symbols binary

 

end n

return

'''Output:'''

<pre>

=={{header|Nim}}==

{{trans|Python}}

 

 

 

that = biased()

 

for n in 3..6:

var cnt0, cnt1 = 0

for x in v:

 

v = newSeqWith(1_000_000, unbiased(biased))

cnt1 = 0

for x in v:

 

=={{header|OCaml}}==

 

if Random.int n = 0 then 1 else 0

 

Printf.printf "%d: %5.2f%% %5.2f%%\n"

b (100.0 *. float !cb /. float n) (100.0 *. float !cu /. float n)

 

Output:

=={{header|PARI/GP}}==

GP's random number generation is high-quality, using Brent's [http://maths.anu.edu.au/~brent/random.html 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 <code>random</code>, 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...

unbiased(N)={

my(a,b);

)

};

 

Output:

 

=={{header|Perl}}==

my $n = shift;

return int(rand($n) / ($n - 1));

100 * sqrt($fixed[0] * $fixed[1]) / ($fixed[0] + $fixed[1])**1.5);

 

Output:

<pre>Bias 3: 6684 3316, 66.84+-0.471% fixed: 2188 2228, 49.5471+-0.752%

Bias 6: 8309 1691, 83.09+-0.375% fixed: 1403 1410, 49.8756+-0.943%</pre>

 

 

=={{header|PicoLisp}}==

(if (= 1 (rand 1 N)) 1 0) )

 

(setq A (randN N))

(setq B (randN N)) ) )

Test:

(tab (2 1 7 2 7 2)

N ":"

(let S 0 (do 10000 (inc 'S (unbiased N))))

2 )

Output:

<pre> 3: 33.21 % 50.48 %

 

=={{header|PL/I}}==

test: procedure options (main); /* 20 Nov. 2012 */

 

 

end test;

Results:

<pre>

5 16 41

6 11 53

</pre>

 

=={{header|PureBasic}}==

If Random(n) <> 1

ProcedureReturn 0

output + #tab$ + " ratio=" + StrF(u_count(1) / #count * 100, 2) + "%" + #LF$

Next

Sample output:

<pre>---------------------------

 

=={{header|Python}}==

import random

 

v = [unbiased(biased) for x in range(1000000)]

v1, v0 = v.count(1), v.count(0)

 

'''Sample output'''

Biased(6) = Stats(count1=167561, count0=832439, percent=16.7561)

Unbiased = Stats(count1=499963, count0=500037, percent=49.996299999999998)</pre>

 

=={{header|R}}==

 

 

unbiased = function(f)

N = N,

biased = mean(replicate(samples, randN(N))),

 

Sample output:

 

=={{header|Racket}}==

#lang racket

;; Using boolean #t/#f instead of 1/0

(printf "Count: ~a => Biased: ~a%; Unbiased: ~a%.\n"

n (try% biased) (try% (unbiased biased))))

{{out}}

<pre>

Count: 6 => Biased: 17%; Unbiased: 50%.

</pre>

 

=={{header|REXX}}==

dash=

end /*N*/

<pre>

</pre>

<pre>

</pre>

<pre>

</pre>

 

=={{header|Ruby}}==

rand(bias) == 0 ? 1 : 0

end

counter[:bias] = bias

puts counter.values_at(*keys).join("\t")

{{output}}

<pre>

5 199767 500354

6 166163 499809

</pre>

 

=={{header|Scala}}==

 

def unbiased( n:Int ) = { def loop : Boolean = { val a = biased(n); if( a != biased(n) ) a else loop }; loop }

"%d: %2.2f%% %2.2f%%".format(i, 100.0*c1/m, 100.0*c2/m)

{{output}}

<pre>3: 33.09% 49.79%

 

=={{header|Seed7}}==

include "float.s7i";

 

" " <& flt(100 * sumUnbiased) / flt(tests) digits 3 lpad 6);

end for;

 

Output:

5: 20.186 49.978

6: 16.570 49.936

</pre>

 

=={{header|Tcl}}==

proc randN n {expr {rand()*$n < 1}}

 

puts [format "unbiased %d => #0=%d #1=%d ratio=%.2f%%" $n $c(0) $c(1) \

[expr {100.*$c(1)/$i}]]

Sample output:

<pre>

biased 6 => #0=833623 #1=166377 ratio=16.64%

unbiased 6 => #0=500518 #1=499482 ratio=49.95%

</pre>

 

=={{header|zkl}}==

foreach N in ([3..6]){

"%d: biased: %3.2f%%, unbiased: %3.2f%%".fmt(N,

(0).reduce(Z,'wrap(s,_){ s+unbiased(randN.fp(N)) },0.0)/Z*100)

.println();

{{out}}

<pre>