Unbias a random generator: Difference between revisions

{{task}}Given a weighted one -bit generator of random numbers where the probability of a one occuringoccurring, <math>P_1</math>, is not the same as <math>P_0</math>, the probability of a zero occuringoccurring, the probability of the occurrence of a one followed by a zero is <math>P_1</math> × <math>P_0</math>, assuming independence. This is the same as the probability of a zero followed by a one: <math>P_0</math> × <math>P_1</math>.

<langsyntaxhighlight Adalang="ada">with Ada.Text_IO; with Ada.Numerics.Discrete_Random;

end Bias_Unbias;</langsyntaxhighlight>

<langsyntaxhighlight AutoHotkeylang="autohotkey">Biased(){

MsgBox % "biased probability of a 1 occurring: " Errorlevel/1000</langsyntaxhighlight>

<langsyntaxhighlight Clang="c">#include <stdio.h>

return rand()r < rand_max / bias;

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

<syntaxhighlight lang c sharp="csharp">using System;

}</langsyntaxhighlight>

<langsyntaxhighlight Clojurelang="clojure">(defn biased [n]

[6 0.91526 0.5])</langsyntaxhighlight>

<langsyntaxhighlight lang="lisp">(defun biased (n) (if (zerop (random n)) 0 1))

(format t "~a: unbiased ~d biased ~d~%" n (count 0 u) (count 0 b))))</langsyntaxhighlight>output<lang>3: unbiased 4992 biased 3361

6: unbiased 4913 biased 1658</langpre>

<langsyntaxhighlight lang="d">import std.stdio, std.random, std.algorithm, std.range, std.functional;

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

boolint unbiased(in int nbias) /*nothrow*/ {

boolint a, b;

intenum M = 50_000500_000;

foreach (immutable n; 3 .. 7) {

intwritefln("%d: c1%2.3f%% %2.3f%%", c2;n,

foreach (i; 0 M.iota.map!(_=> Mn.biased).sum {* 100.0 / M,

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

<langsyntaxhighlight lang="euphoria">function randN(integer N)

end for</langsyntaxhighlight>

<langsyntaxhighlight lang="fortran">program Bias_Unbias

end program</langsyntaxhighlight>

<langsyntaxhighlight lang="gap">RandNGen := function(n)

# [ 6, 166746, 499859 ] ]</langsyntaxhighlight>

<langsyntaxhighlight lang="go">package main

"math/rand"

}</langsyntaxhighlight>

<langsyntaxhighlight lang="haskell">import Control.Monad.Random

import RandomControl.Monad

randN :: IntegerMonadRandom m => Int -> IOm BoolInt

randN n = randomRIOfromList [(10, fromIntegral n-1) >>= return ., (==1, 1)]</syntaxhighlight>

λ> replicateM b <-20 (randN n5)

where trials = 50000</lang>

n = 3 biased: 33.72% 50.08unbiased: 49%

n = 4 biased: 25.2624% unbiased: 50.15%

n = 5 biased: 19.99% unbiased: 50.07%

n = 6 biased: 16.67% 50.10unbiased: 49%

<langsyntaxhighlight lang="j">randN=: 0 = ?

unbiased=: i.@# { ::$: 2 | 0 3 -.~ _2 #.\ 4&* randN@# ]</langsyntaxhighlight>

<langsyntaxhighlight lang="j"> randN 10#6

1 0 0 1 0 0 1 0 1 1</langsyntaxhighlight>

<langsyntaxhighlight lang="j"> +/randN 100#3

47</langsyntaxhighlight>

Note that these results are randomarbitrary. For example, a re-run of <code>+/randN 100#5</code> gave 25 as its result, and a re-run of <code>+/unbiased 100#5</code> gave 52 as its result.

<langsyntaxhighlight lang="java">public class Bias {

}</langsyntaxhighlight>

<langsyntaxhighlight lang="lua">

<langsyntaxhighlight lang="ocaml">let randN n =

done</langsyntaxhighlight>

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...

<langsyntaxhighlight lang="parigp">randN(N)=!random(N);

for(n=3,6,print(n"\t"sum(k=1,1e5,unbiased(n))"\t"sum(k=1,1e5,randN(n))))</langsyntaxhighlight>

<langsyntaxhighlight lang="perl">sub randn {

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

=={{header|Perl 6Phix}}==

<langsyntaxhighlight PicoLisplang="picolisp">(de randN (N)

A ) )</langsyntaxhighlight>

<langsyntaxhighlight PicoLisplang="picolisp">(for N (range 3 6)

"%" ) )</langsyntaxhighlight>

<langsyntaxhighlight PureBasiclang="purebasic">Procedure biased(n)

MessageRequester("Biased and Unbiased random number results", output)</langsyntaxhighlight>

<langsyntaxhighlight lang="python">from __future__ import print_function

print ( " Unbiased = %r" % (Stats(v1, v0, 100. * v1/(v1 + v0)), ) )</langsyntaxhighlight>

<langsyntaxhighlight lang="seed7">$ include "seed7_05.s7i";

end func;</langsyntaxhighlight>

<langsyntaxhighlight lang="tcl"># 1,0 random generator factory with 1 appearing 1/N'th of the time

}</langsyntaxhighlight>