Rosetta Code

Percolation/Mean run density: Difference between revisions

Page Discussion

Percolation/Mean run density (view source)

Revision as of 11:58, 24 January 2024

30,345 bytes added ,  3 months ago
m
→‎{{header|Wren}}: Minor tidy
(Added EchoLisp)
m (→‎{{header|Wren}}: Minor tidy)
 
(29 intermediate revisions by 17 users not shown)
Line 27:
;See also
* [http://mathworld.wolfram.com/s-Run.html s-Run] on Wolfram mathworld.
 
=={{header|11l}}==
{{trans|Python}}
 
<syntaxhighlight lang="11l">UInt32 seed = 0
F nonrandom()
:seed = 1664525 * :seed + 1013904223
R Int(:seed >> 16) / Float(FF'FF)
 
V (p, t) = (0.5, 500)
 
F newv(n, p)
R (0 .< n).map(i -> Int(nonrandom() < @p))
 
F runs(v)
R sum(zip(v, v[1..] [+] [0]).map((a, b) -> (a [&] ~b)))
 
F mean_run_density(n, p)
R runs(newv(n, p)) / Float(n)
 
L(p10) (1.<10).step(2)
p = p10 / 10
V limit = p * (1 - p)
print(‘’)
L(n2) (10.<16).step(2)
V n = 2 ^ n2
V sim = sum((0 .< t).map(i -> mean_run_density(@n, :p))) / t
print(‘t=#3 p=#.2 n=#5 p(1-p)=#.3 sim=#.3 delta=#.1%’.format(
t, p, n, limit, sim, I limit {abs(sim - limit) / limit * 100} E sim * 100))</syntaxhighlight>
 
{{out}}
<pre>
 
t=500 p=0.10 n= 1024 p(1-p)=0.090 sim=0.090 delta=0.0%
t=500 p=0.10 n= 4096 p(1-p)=0.090 sim=0.090 delta=0.1%
t=500 p=0.10 n=16384 p(1-p)=0.090 sim=0.090 delta=0.0%
 
t=500 p=0.30 n= 1024 p(1-p)=0.210 sim=0.210 delta=0.1%
t=500 p=0.30 n= 4096 p(1-p)=0.210 sim=0.210 delta=0.0%
t=500 p=0.30 n=16384 p(1-p)=0.210 sim=0.210 delta=0.0%
 
t=500 p=0.50 n= 1024 p(1-p)=0.250 sim=0.251 delta=0.2%
t=500 p=0.50 n= 4096 p(1-p)=0.250 sim=0.250 delta=0.1%
t=500 p=0.50 n=16384 p(1-p)=0.250 sim=0.250 delta=0.0%
 
t=500 p=0.70 n= 1024 p(1-p)=0.210 sim=0.211 delta=0.3%
t=500 p=0.70 n= 4096 p(1-p)=0.210 sim=0.210 delta=0.1%
t=500 p=0.70 n=16384 p(1-p)=0.210 sim=0.210 delta=0.0%
 
t=500 p=0.90 n= 1024 p(1-p)=0.090 sim=0.091 delta=1.0%
t=500 p=0.90 n= 4096 p(1-p)=0.090 sim=0.090 delta=0.0%
t=500 p=0.90 n=16384 p(1-p)=0.090 sim=0.090 delta=0.0%
</pre>
 
=={{header|ALGOL 68}}==
{{Trans|C}}
<syntaxhighlight lang="algol68">
BEGIN
 
# just generate 0s and 1s without storing them #
PROC run test = ( REAL p, INT len, runs )REAL:
BEGIN
INT count := 0;
REAL thresh = p;
TO runs DO
INT x := 0;
FOR i FROM len BY -1 TO 1 DO
INT y = ABS ( random < thresh );
count +:= ABS ( x < y );
x := y
OD
OD;
count / runs / len
END # run test # ;
 
print( ( "running 1000 tests each:", newline ) );
print( ( " p n K p(1-p) diff", newline ) );
print( ( "----------------------------------------------", newline ) );
FOR ip BY 2 TO 9 DO
REAL p = ip / 10;
REAL p1p = p * (1 - p);
INT n := 10;
WHILE ( n *:= 10 ) <= 100000 DO
REAL k = run test( p, n, 1000 );
print( ( fixed( p, -4, 1 ), whole( n, -9 ), fixed( k, -8, 4 )
, fixed( p1p, -8, 4 ), fixed( k - p1p, 9, 4 )
, " (", fixed( ( k - p1p ) / p1p * 100, 5, 2 ), "%)", newline
)
)
OD;
print( ( newline ) )
OD
END
</syntaxhighlight>
{{out}}
<pre>
running 1000 tests each:
p n K p(1-p) diff
----------------------------------------------
0.1 100 0.0898 0.0900 -0.0002 (-0.21%)
0.1 1000 0.0902 0.0900 +0.0002 (+0.20%)
0.1 10000 0.0900 0.0900 +0.0000 (+0.05%)
0.1 100000 0.0900 0.0900 +0.0000 (+0.04%)
 
0.3 100 0.2105 0.2100 +0.0005 (+0.22%)
0.3 1000 0.2098 0.2100 -0.0002 (-0.12%)
0.3 10000 0.2100 0.2100 +0.0000 (+0.02%)
0.3 100000 0.2101 0.2100 +0.0001 (+0.03%)
 
0.5 100 0.2536 0.2500 +0.0035 (+1.42%)
0.5 1000 0.2504 0.2500 +0.0004 (+0.16%)
0.5 10000 0.2501 0.2500 +0.0001 (+0.03%)
0.5 100000 0.2500 0.2500 +0.0000 (+0.01%)
 
0.7 100 0.2155 0.2100 +0.0055 (+2.60%)
0.7 1000 0.2107 0.2100 +0.0007 (+0.33%)
0.7 10000 0.2101 0.2100 +0.0001 (+0.06%)
0.7 100000 0.2100 0.2100 -0.0000 (-0.02%)
 
0.9 100 0.0982 0.0900 +0.0082 (+9.12%)
0.9 1000 0.0902 0.0900 +0.0002 (+0.27%)
0.9 10000 0.0901 0.0900 +0.0001 (+0.11%)
0.9 100000 0.0900 0.0900 +0.0000 (+0.01%)
</pre>
 
=={{header|C}}==
<langsyntaxhighlight lang="c">#include <stdio.h>
#include <stdlib.h>
 
Line 64 ⟶ 188:
 
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>
Line 95 ⟶ 219:
0.9 100000 0.0900 0.0900 +0.0000 (+0.03%)
</pre>
 
=={{header|C++}}==
<langsyntaxhighlight lang="cpp">#include <algorithm>
#include <random>
#include <vector>
Line 151 ⟶ 276:
}
return 0 ;
}</langsyntaxhighlight>
{{out}}
<pre>t = 100
Line 183 ⟶ 308:
=={{header|D}}==
{{trans|python}}
<langsyntaxhighlight lang="d">import std.stdio, std.range, std.algorithm, std.random, std.math;
 
enum n = 100, p = 0.5, t = 500;
Line 205 ⟶ 330:
}
}
}</langsyntaxhighlight>
{{out}}
<pre>t=500, p=0.10, n= 1024, p(1-p)=0.09000, sim=0.08949, delta=0.6%
Line 226 ⟶ 351:
t=500, p=0.90, n= 4096, p(1-p)=0.09000, sim=0.09047, delta=0.5%
t=500, p=0.90, n=16384, p(1-p)=0.09000, sim=0.09007, delta=0.1%</pre>
 
=={{header|EasyLang}}==
<syntaxhighlight lang="easylang">
numfmt 3 6
for p in [ 0.1 0.3 0.5 0.7 0.9 ]
theory = p * (1 - p)
print "p:" & p & " theory:" & theory
print " n sim"
for n in [ 1e2 1e3 1e4 ]
sum = 0
for t to 100
run = 0
for j to n
h = if randomf < p
if h = 1 and run = 0
sum += 1
.
run = h
.
.
print n & " " & sum / n / t
.
print ""
.
</syntaxhighlight>
 
=={{header|EchoLisp}}==
<langsyntaxhighlight lang="scheme">
;; count 1-runs - The vector is not stored
(define (runs p n)
Line 251 ⟶ 401:
(for ([n '(10 100 1000)])
(printf "\t-- n %5d → %d" n (truns p n)))))
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 275 ⟶ 425:
-- n 100 → 0.0894
-- n 1000 → 0.0905
</pre>
 
=={{header|Factor}}==
<syntaxhighlight lang="factor">USING: formatting fry io kernel math math.ranges math.statistics
random sequences ;
IN: rosetta-code.mean-run-density
 
: rising? ( ? ? -- ? ) [ f = ] [ t = ] bi* and ;
 
: count-run ( n ? ? -- m ? )
2dup rising? [ [ 1 + ] 2dip ] when nip ;
 
: runs ( n p -- n )
[ 0 f ] 2dip '[ random-unit _ < count-run ] times drop ;
 
: rn ( n p -- x ) over [ runs ] dip /f ;
 
: sim ( n p -- avg )
[ 1000 ] 2dip [ rn ] 2curry replicate mean ;
 
: theory ( p -- x ) 1 over - * ;
 
: result ( n p -- )
[ swap ] [ sim ] [ nip theory ] 2tri 2dup - abs
"%.1f %-5d %.4f %.4f %.4f\n" printf ;
 
: test ( p -- )
{ 100 1,000 10,000 } [ swap result ] with each nl ;
 
: header ( -- )
"1000 tests each:\np n K p(1-p) diff" print ;
 
: main ( -- ) header .1 .9 .2 <range> [ test ] each ;
 
MAIN: main</syntaxhighlight>
{{out}}
<pre>
1000 tests each:
p n K p(1-p) diff
0.1 100 0.0909 0.0900 0.0009
0.1 1000 0.0902 0.0900 0.0002
0.1 10000 0.0899 0.0900 0.0001
 
0.3 100 0.2111 0.2100 0.0011
0.3 1000 0.2101 0.2100 0.0001
0.3 10000 0.2100 0.2100 0.0000
 
0.5 100 0.2524 0.2500 0.0024
0.5 1000 0.2504 0.2500 0.0004
0.5 10000 0.2501 0.2500 0.0001
 
0.7 100 0.2149 0.2100 0.0049
0.7 1000 0.2106 0.2100 0.0006
0.7 10000 0.2100 0.2100 0.0000
 
0.9 100 0.0978 0.0900 0.0078
0.9 1000 0.0905 0.0900 0.0005
0.9 10000 0.0901 0.0900 0.0001
</pre>
 
=={{header|Fortran}}==
<langsyntaxhighlight lang="fortran">
! loosely translated from python. We do not need to generate and store the entire vector at once.
! compilation: gfortran -Wall -std=f2008 -o thisfile thisfile.f08
Line 338 ⟶ 546:
 
end program percolation_mean_run_density
</syntaxhighlight>
</lang>
 
<pre>
Line 364 ⟶ 572:
500 0.90 16384 0.090 0.090 0.1
</pre>
 
=={{header|FreeBASIC}}==
{{trans|Phix}}
<syntaxhighlight lang="freebasic">Function run_test(p As Double, longitud As Integer, runs As Integer) As Double
Dim As Integer r, l, cont = 0
Dim As Integer v, pv
For r = 1 To runs
pv = 0
For l = 1 To longitud
v = Rnd < p
cont += Iif(pv < v, 1, 0)
pv = v
Next l
Next r
Return (cont/runs/longitud)
End Function
 
Print "Running 1000 tests each:"
Print " p n K p(1-p) delta"
Print String(46,"-")
 
Dim As Double K, p, p1p
Dim As Integer n, ip
 
For ip = 1 To 10 Step 2
p = ip / 10
p1p = p * (1-p)
n = 100
While n <= 100000
K = run_test(p, n, 1000)
Print Using !"#.# ###### #.#### #.#### +##.#### (##.## \b%)"; _
p; n; K; p1p; K-p1p; (K-p1p)/p1p*100
n *= 10
Wend
Print
Next ip
Sleep</syntaxhighlight>
<pre>Same as Phix, C, Kotlin, Wren, Pascal or zkl entry.</pre>
 
=={{header|Go}}==
<langsyntaxhighlight lang="go">package main
 
import (
Line 400 ⟶ 647:
}
}
}</langsyntaxhighlight>
{{out}}
<pre>
Line 440 ⟶ 687:
 
=={{header|Haskell}}==
<langsyntaxhighlight Haskelllang="haskell">import Control.Monad.Random
import Control.Applicative
import Text.Printf
Line 468 ⟶ 715:
x <- newStdGen
forM_ [0.1,0.3,0.5,0.7,0.9] $ printKs x
</syntaxhighlight>
</lang>
<pre>./percolation
p= 0.1, K(p)= 0.090
Line 506 ⟶ 753:
The following works in both languages:
 
<langsyntaxhighlight lang="unicon">procedure main(A)
t := integer(A[2]) | 500
 
Line 520 ⟶ 767:
write(left(p,8)," ",left(n,8)," ",left(p*(1-p),10)," ",left(Ka/t, 10))
}
end</langsyntaxhighlight>
 
Output:
Line 546 ⟶ 793:
 
=={{header|J}}==
<syntaxhighlight lang="j">
<lang J>
NB. translation of python
NB. 'N P T' =: 100 0.5 500 NB. sillinesshypothetical example values, to aid comprehension...
newv =: (> ?@(#&0))~ NB. generate a random binary vector. Use: N newv P
Line 568 ⟶ 815:
EMPTY
)
</syntaxhighlight>
</lang>
Session:
<pre>
Line 595 ⟶ 842:
</pre>
 
=={{header|MathematicaJava}}==
<syntaxhighlight lang="java">
<lang Mathematica>meanRunDensity[p_, len_, trials_] :=
 
import java.util.concurrent.ThreadLocalRandom;
 
public final class PercolationMeanRun {
 
public static void main(String[] aArgs) {
System.out.println("Running 1000 tests each:" + System.lineSeparator());
System.out.println(" p\tlength\tresult\ttheory\t difference");
System.out.println("-".repeat(48));
for ( double probability = 0.1; probability <= 0.9; probability += 0.2 ) {
double theory = probability * ( 1.0 - probability );
int length = 100;
while ( length <= 100_000 ) {
double result = runTest(probability, length, 1_000);
System.out.println(String.format("%.1f\t%6d\t%.4f\t%.4f\t%+.4f (%+.2f%%)",
probability, length, result, theory, result - theory, ( result - theory ) / theory * 100));
length *= 10;
}
System.out.println();
}
 
}
private static double runTest(double aProbability, int aLength, int aRunCount) {
double count = 0.0;
for ( int run = 0; run < aRunCount; run++ ) {
int previousBit = 0;
int length = aLength;
while ( length-- > 0 ) {
int nextBit = ( random.nextDouble(1.0) < aProbability ) ? 1 : 0;
if ( previousBit < nextBit ) {
count += 1.0;
}
previousBit = nextBit;
}
}
return count / aRunCount / aLength;
}
 
private static ThreadLocalRandom random = ThreadLocalRandom.current();
 
}
</syntaxhighlight>
{{ out }}
<pre>
Running 1000 tests each:
 
p length result theory difference
------------------------------------------------
0.1 100 0.0899 0.0900 -0.0001 (-0.07%)
0.1 1000 0.0902 0.0900 +0.0002 (+0.18%)
0.1 10000 0.0900 0.0900 +0.0000 (+0.02%)
0.1 100000 0.0900 0.0900 -0.0000 (-0.00%)
 
0.3 100 0.2110 0.2100 +0.0010 (+0.47%)
0.3 1000 0.2101 0.2100 +0.0001 (+0.05%)
0.3 10000 0.2100 0.2100 -0.0000 (-0.01%)
0.3 100000 0.2100 0.2100 -0.0000 (-0.01%)
 
0.5 100 0.2516 0.2500 +0.0015 (+0.62%)
0.5 1000 0.2509 0.2500 +0.0009 (+0.37%)
0.5 10000 0.2499 0.2500 -0.0001 (-0.04%)
0.5 100000 0.2500 0.2500 +0.0000 (+0.00%)
 
0.7 100 0.2145 0.2100 +0.0045 (+2.12%)
0.7 1000 0.2106 0.2100 +0.0006 (+0.28%)
0.7 10000 0.2101 0.2100 +0.0001 (+0.06%)
0.7 100000 0.2100 0.2100 -0.0000 (-0.00%)
 
0.9 100 0.0970 0.0900 +0.0070 (+7.74%)
0.9 1000 0.0910 0.0900 +0.0010 (+1.15%)
0.9 10000 0.0901 0.0900 +0.0001 (+0.06%)
0.9 100000 0.0900 0.0900 +0.0000 (+0.00%)
</pre>
 
=={{header|Julia}}==
{{trans|Python}}
<syntaxhighlight lang="julia">using Printf, Distributions, IterTools
 
newv(n::Int, p::Float64) = rand(Bernoulli(p), n)
runs(v::Vector{Int}) = sum((a & ~b) for (a, b) in zip(v, IterTools.chain(v[2:end], v[1])))
 
mrd(n::Int, p::Float64) = runs(newv(n, p)) / n
 
nrep = 500
 
for p in 0.1:0.2:1
lim = p * (1 - p)
 
println()
for ex in 10:2:14
n = 2 ^ ex
sim = mean(mrd.(n, p) for _ in 1:nrep)
@printf("nrep = %3i\tp = %4.2f\tn = %5i\np · (1 - p) = %5.3f\tsim = %5.3f\tΔ = %3.1f%%\n",
nrep, p, n, lim, sim, lim > 0 ? abs(sim - lim) / lim * 100 : sim * 100)
end
end</syntaxhighlight>
 
{{out}}
<pre>
nrep = 500 p = 0.10 n = 1024
p · (1 - p) = 0.090 sim = 0.090 Δ = 0.4%
nrep = 500 p = 0.10 n = 4096
p · (1 - p) = 0.090 sim = 0.090 Δ = 0.2%
nrep = 500 p = 0.10 n = 16384
p · (1 - p) = 0.090 sim = 0.090 Δ = 0.0%
 
nrep = 500 p = 0.30 n = 1024
p · (1 - p) = 0.210 sim = 0.211 Δ = 0.5%
nrep = 500 p = 0.30 n = 4096
p · (1 - p) = 0.210 sim = 0.210 Δ = 0.1%
nrep = 500 p = 0.30 n = 16384
p · (1 - p) = 0.210 sim = 0.210 Δ = 0.0%
 
nrep = 500 p = 0.50 n = 1024
p · (1 - p) = 0.250 sim = 0.250 Δ = 0.0%
nrep = 500 p = 0.50 n = 4096
p · (1 - p) = 0.250 sim = 0.250 Δ = 0.1%
nrep = 500 p = 0.50 n = 16384
p · (1 - p) = 0.250 sim = 0.250 Δ = 0.0%
 
nrep = 500 p = 0.70 n = 1024
p · (1 - p) = 0.210 sim = 0.209 Δ = 0.3%
nrep = 500 p = 0.70 n = 4096
p · (1 - p) = 0.210 sim = 0.210 Δ = 0.1%
nrep = 500 p = 0.70 n = 16384
p · (1 - p) = 0.210 sim = 0.210 Δ = 0.0%
 
nrep = 500 p = 0.90 n = 1024
p · (1 - p) = 0.090 sim = 0.090 Δ = 0.0%
nrep = 500 p = 0.90 n = 4096
p · (1 - p) = 0.090 sim = 0.090 Δ = 0.0%
nrep = 500 p = 0.90 n = 16384
p · (1 - p) = 0.090 sim = 0.090 Δ = 0.1%</pre>
 
=={{header|Kotlin}}==
{{trans|C}}
<syntaxhighlight lang="scala">// version 1.2.10
 
import java.util.Random
 
val rand = Random()
const val RAND_MAX = 32767
 
// just generate 0s and 1s without storing them
fun runTest(p: Double, len: Int, runs: Int): Double {
var cnt = 0
val thresh = (p * RAND_MAX).toInt()
for (r in 0 until runs) {
var x = 0
var i = len
while (i-- > 0) {
val y = if (rand.nextInt(RAND_MAX + 1) < thresh) 1 else 0
if (x < y) cnt++
x = y
}
}
return cnt.toDouble() / runs / len
}
 
fun main(args: Array<String>) {
println("running 1000 tests each:")
println(" p\t n\tK\tp(1-p)\t diff")
println("------------------------------------------------")
val fmt = "%.1f\t%6d\t%.4f\t%.4f\t%+.4f (%+.2f%%)"
for (ip in 1..9 step 2) {
val p = ip / 10.0
val p1p = p * (1.0 - p)
var n = 100
while (n <= 100_000) {
val k = runTest(p, n, 1000)
println(fmt.format(p, n, k, p1p, k - p1p, (k - p1p) / p1p * 100))
n *= 10
}
println()
}
}</syntaxhighlight>
 
Sample output:
<pre>
running 1000 tests each:
p n K p(1-p) diff
------------------------------------------------
0.1 100 0.0908 0.0900 +0.0008 (+0.93%)
0.1 1000 0.0900 0.0900 +0.0000 (+0.02%)
0.1 10000 0.0899 0.0900 -0.0001 (-0.08%)
0.1 100000 0.0900 0.0900 -0.0000 (-0.05%)
 
0.3 100 0.2112 0.2100 +0.0012 (+0.56%)
0.3 1000 0.2096 0.2100 -0.0004 (-0.21%)
0.3 10000 0.2101 0.2100 +0.0001 (+0.05%)
0.3 100000 0.2101 0.2100 +0.0001 (+0.03%)
 
0.5 100 0.2522 0.2500 +0.0022 (+0.90%)
0.5 1000 0.2504 0.2500 +0.0004 (+0.15%)
0.5 10000 0.2500 0.2500 -0.0000 (-0.00%)
0.5 100000 0.2500 0.2500 +0.0000 (+0.00%)
 
0.7 100 0.2162 0.2100 +0.0062 (+2.95%)
0.7 1000 0.2106 0.2100 +0.0006 (+0.29%)
0.7 10000 0.2101 0.2100 +0.0001 (+0.03%)
0.7 100000 0.2100 0.2100 +0.0000 (+0.01%)
 
0.9 100 0.0982 0.0900 +0.0083 (+9.17%)
0.9 1000 0.0911 0.0900 +0.0011 (+1.17%)
0.9 10000 0.0902 0.0900 +0.0002 (+0.18%)
0.9 100000 0.0900 0.0900 -0.0000 (-0.02%)
</pre>
 
=={{header|Mathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">meanRunDensity[p_, len_, trials_] :=
Mean[Length[Cases[Split@#, {1, ___}]] & /@
Unitize[Chop[RandomReal[1, {trials, len}], 1 - p]]]/len
Line 604 ⟶ 1,063:
Table[{q, n, x = meanRunDensity[q, n, 100] // N,
q (1 - q) - x}, {n, {100, 1000, 10000, 100000}}], {}],
Alignment -> Left], {q, {.1, .3, .5, .7, .9}}]</langsyntaxhighlight>
{{out}}
<pre>
Line 637 ⟶ 1,096:
0.9 100000 0.0900 -0.0000144
</pre>
 
=={{header|Nim}}==
{{trans|Go}}
<syntaxhighlight lang="nim">import random, strformat
 
const T = 100
 
var
pList = [0.1, 0.3, 0.5, 0.7, 0.9]
nList = [100, 1_000, 10_000, 100_000]
 
for p in pList:
 
let theory = p * (1 - p)
echo &"\np: {p:.4f} theory: {theory:.4f} t: {T}"
echo " n sim sim-theory"
 
for n in nList:
var sum = 0
for _ in 1..T:
var run = false
for _ in 1..n:
let one = rand(1.0) < p
if one and not run: inc sum
run = one
 
let k = sum / (T * n)
echo &"{n:9} {k:15.4f} {k - theory:10.6f}"</syntaxhighlight>
 
{{out}}
<pre>p: 0.1000 theory: 0.0900 t: 100
n sim sim-theory
100 0.0886 -0.001400
1000 0.0907 0.000750
10000 0.0903 0.000330
100000 0.0900 -0.000013
 
p: 0.3000 theory: 0.2100 t: 100
n sim sim-theory
100 0.2135 0.003500
1000 0.2086 -0.001390
10000 0.2102 0.000180
100000 0.2099 -0.000132
 
p: 0.5000 theory: 0.2500 t: 100
n sim sim-theory
100 0.2546 0.004600
1000 0.2495 -0.000550
10000 0.2502 0.000232
100000 0.2500 -0.000032
 
p: 0.7000 theory: 0.2100 t: 100
n sim sim-theory
100 0.2148 0.004800
1000 0.2110 0.001010
10000 0.2105 0.000525
100000 0.2099 -0.000077
 
p: 0.9000 theory: 0.0900 t: 100
n sim sim-theory
100 0.0968 0.006800
1000 0.0916 0.001570
10000 0.0903 0.000334
100000 0.0901 0.000113</pre>
 
=={{header|Pascal}}==
{{trans|C}}{{works with|Free Pascal}}
<syntaxhighlight lang="pascal">
{$MODE objFPC}//for using result,parameter runs becomes for variable..
uses
sysutils;//Format
const
MaxN = 100*1000;
 
function run_test(p:double;len,runs: NativeInt):double;
var
x, y, i,cnt : NativeInt;
Begin
result := 1/ (runs * len);
cnt := 0;
for runs := runs-1 downto 0 do
Begin
x := 0;
y := 0;
for i := len-1 downto 0 do
begin
x := y;
y := Ord(Random() < p);
cnt := cnt+ord(x < y);
end;
end;
result := result *cnt;
end;
 
//main
var
p, p1p, K : double;
ip, n : nativeInt;
Begin
randomize;
writeln( 'running 1000 tests each:'#13#10,
' p n K p(1-p) diff'#13#10,
'-----------------------------------------------');
ip:= 1;
while ip < 10 do
Begin
p := ip / 10;
p1p := p * (1 - p);
n := 100;
While n <= MaxN do
Begin
K := run_test(p, n, 1000);
writeln(Format('%4.1f %6d %6.4f %6.4f %7.4f (%5.2f %%)',
[p, n, K, p1p, K - p1p, (K - p1p) / p1p * 100]));
n := n*10;
end;
writeln;
ip := ip+2;
end;
end.</syntaxhighlight>
Output
<pre>running 1000 tests each:
p n K p(1-p) diff
-----------------------------------------------
0.1 100 0.0894 0.0900 -0.0006 (-0.70 %)
0.1 1000 0.0898 0.0900 -0.0002 (-0.17 %)
0.1 10000 0.0900 0.0900 0.0000 ( 0.02 %)
0.1 100000 0.0900 0.0900 0.0000 ( 0.04 %)
 
0.3 100 0.2112 0.2100 0.0012 ( 0.57 %)
0.3 1000 0.2101 0.2100 0.0001 ( 0.04 %)
0.3 10000 0.2099 0.2100 -0.0001 (-0.04 %)
0.3 100000 0.2099 0.2100 -0.0001 (-0.03 %)
 
0.5 100 0.2516 0.2500 0.0016 ( 0.66 %)
0.5 1000 0.2497 0.2500 -0.0003 (-0.14 %)
0.5 10000 0.2501 0.2500 0.0001 ( 0.03 %)
0.5 100000 0.2500 0.2500 0.0000 ( 0.01 %)
 
0.7 100 0.2144 0.2100 0.0044 ( 2.08 %)
0.7 1000 0.2107 0.2100 0.0007 ( 0.32 %)
0.7 10000 0.2101 0.2100 0.0001 ( 0.02 %)
0.7 100000 0.2100 0.2100 0.0000 ( 0.01 %)
 
0.9 100 0.0978 0.0900 0.0078 ( 8.69 %)
0.9 1000 0.0909 0.0900 0.0009 ( 0.96 %)
0.9 10000 0.0901 0.0900 0.0001 ( 0.10 %)
0.9 100000 0.0900 0.0900 0.0000 ( 0.02 %)
</pre>
 
=={{header|Perl}}==
{{trans|Perl 6Raku}}
<langsyntaxhighlight lang="perl">sub R {
my ($n, $p) = @_;
my $r = join '',
Line 655 ⟶ 1,264:
printf " R(n, p)= %f\n", $r / t;
}
}</langsyntaxhighlight>
{{out}}
<pre>t= 100
Line 679 ⟶ 1,288:
R(n, p)= 0.091730</pre>
 
=={{header|Perl 6Phix}}==
{{trans|zkl}}
<lang perl6>sub R($n, $p) { [+] ((rand < $p) xx $n).squish }
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">run_test</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">integer</span> <span style="color: #000000;">len</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">runs</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">count</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">r</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">runs</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">bool</span> <span style="color: #000000;">v</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">pv</span> <span style="color: #0000FF;">=</span> <span style="color: #004600;">false</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">l</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">len</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">v</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">rnd</span><span style="color: #0000FF;">()<</span><span style="color: #000000;">p</span>
<span style="color: #000000;">count</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">pv</span><span style="color: #0000FF;"><</span><span style="color: #000000;">v</span>
<span style="color: #000000;">pv</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">v</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">return</span> <span style="color: #000000;">count</span><span style="color: #0000FF;">/</span><span style="color: #000000;">runs</span><span style="color: #0000FF;">/</span><span style="color: #000000;">len</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Running 1000 tests each:\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" p n K p(1-p) delta\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"--------------------------------------------\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">ip</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">10</span> <span style="color: #008080;">by</span> <span style="color: #000000;">2</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ip</span><span style="color: #0000FF;">/</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span>
<span style="color: #000000;">p1p</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">*(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">-</span><span style="color: #000000;">p</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">100</span>
<span style="color: #008080;">while</span> <span style="color: #000000;">n</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">100000</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">K</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">run_test</span><span style="color: #0000FF;">(</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1000</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%.1f %6d %6.4f %6.4f %+7.4f (%+5.2f%%)\n"</span><span style="color: #0000FF;">,</span>
<span style="color: #0000FF;">{</span><span style="color: #000000;">p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">K</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">p1p</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">K</span><span style="color: #0000FF;">-</span><span style="color: #000000;">p1p</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">K</span><span style="color: #0000FF;">-</span><span style="color: #000000;">p1p</span><span style="color: #0000FF;">)/</span><span style="color: #000000;">p1p</span><span style="color: #0000FF;">*</span><span style="color: #000000;">100</span><span style="color: #0000FF;">})</span>
<span style="color: #000000;">n</span> <span style="color: #0000FF;">*=</span> <span style="color: #000000;">10</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">while</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #000000;">main</span><span style="color: #0000FF;">()</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Running 1000 tests each:
p n K p(1-p) delta
--------------------------------------------
0.1 100 0.0889 0.0900 -0.0011 (-1.20%)
0.1 1000 0.0896 0.0900 -0.0004 (-0.45%)
0.1 10000 0.0900 0.0900 -0.0000 (-0.02%)
0.1 100000 0.0900 0.0900 -0.0000 (-0.04%)
 
0.3 100 0.2112 0.2100 +0.0012 (+0.57%)
say 't= ', constant t = 100;
0.3 1000 0.2101 0.2100 +0.0001 (+0.07%)
0.3 10000 0.2101 0.2100 +0.0001 (+0.06%)
0.3 100000 0.2100 0.2100 -0.0000 (-0.01%)
 
0.5 100 0.2528 0.2500 +0.0028 (+1.13%)
for .1, .3 ... .9 -> $p {
0.5 1000 0.2500 0.2500 +0.0000 (+0.01%)
say "p= $p, K(p)= {$p*(1-$p)}";
0.5 10000 0.2500 0.2500 +0.0000 (+0.00%)
for 10, 100, 1000 -> $n {
0.5 100000 0.2500 0.2500 -0.0000 (-0.00%)
printf " R(%6d, p)= %f\n", $n, t R/ [+] R($n, $p)/$n xx t
 
}
0.7 100 0.2174 0.2100 +0.0074 (+3.50%)
}</lang>
0.7 1000 0.2105 0.2100 +0.0005 (+0.26%)
{{out}}
0.7 10000 0.2101 0.2100 +0.0001 (+0.06%)
<pre>t= 100
0.7 100000 0.2100 0.2100 +0.0000 (+0.01%)
p= 0.1, K(p)= 0.09
 
R( 10, p)= 0.088000
0.9 R( 100, p)= 0.0856000986 0.0900 +0.0086 (+9.53%)
0.9 R( 1000, p)= 0.0891500908 0.0900 +0.0008 (+0.88%)
0.9 10000 0.0901 0.0900 +0.0001 (+0.11%)
p= 0.3, K(p)= 0.21
0.9 R(100000 0.0900 10,0.0900 p)= +0.2110000000 (+0.03%)
</pre>
R( 100, p)= 0.214600
R( 1000, p)= 0.211160
p= 0.5, K(p)= 0.25
R( 10, p)= 0.279000
R( 100, p)= 0.249200
R( 1000, p)= 0.250870
p= 0.7, K(p)= 0.21
R( 10, p)= 0.258000
R( 100, p)= 0.215400
R( 1000, p)= 0.209560
p= 0.9, K(p)= 0.09
R( 10, p)= 0.181000
R( 100, p)= 0.094500
R( 1000, p)= 0.091330</pre>
 
=={{header|Python}}==
<langsyntaxhighlight lang="python">from __future__ import division
from random import random
from math import fsum
Line 737 ⟶ 1,379:
sim = fsum(mean_run_density(n, p) for i in range(t)) / t
print('t=%3i p=%4.2f n=%5i p(1-p)=%5.3f sim=%5.3f delta=%3.1f%%'
% (t, p, n, limit, sim, abs(sim - limit) / limit * 100 if limit else sim * 100))</langsyntaxhighlight>
 
{{out}}
Line 762 ⟶ 1,404:
=={{header|Racket}}==
 
<langsyntaxhighlight lang="racket">#lang racket
(require racket/fixnum)
(define t (make-parameter 100))
Line 793 ⟶ 1,435:
(module+ test
(require rackunit)
(check-eq? (Rn (fxvector 1 1 0 0 0 1 0 1 1 1)) 3))</langsyntaxhighlight>
{{out}}
<pre>
Line 827 ⟶ 1,469:
 
</pre>
 
=={{header|REXX}}==
{{trans|Fortran}}
<syntaxhighlight lang="rexx">/* REXX */
Numeric Digits 20
Call random(,12345) /* make the run reproducable */
pList = '.1 .3 .5 .7 .9'
nList = '1e2 1e3 1e4 1e5'
t = 100
Do While plist<>''
Parse Var plist p plist
theory=p*(1-p)
Say ' '
Say 'p:' format(p,2,4)' theory:'format(theory,2,4)' t:'format(t,4)
Say ' n sim sim-theory'
nl=nlist
Do While nl<>''
Parse Var nl n nl
sum=0
Do i=1 To t
run=0
Do j=1 To n
one=random(1000)<p*1000
If one & (run=0) Then
sum=sum+1
run=one
End
End
sim=sum/(n*100)
Say format(n,10)' ' format(sim,2,4)' 'format(sim-theory,2,6)
End
End</syntaxhighlight>
{{out}}
<pre>p: 0.1000 theory: 0.0900 t: 100
n sim sim-theory
100 0.0875 -0.002500
1000 0.0894 -0.000560
10000 0.0902 0.000237
100000 0.0899 -0.000112
 
p: 0.3000 theory: 0.2100 t: 100
n sim sim-theory
100 0.2088 -0.001200
1000 0.2116 0.001570
10000 0.2101 0.000056
100000 0.2099 -0.000120
 
p: 0.5000 theory: 0.2500 t: 100
n sim sim-theory
100 0.2557 0.005700
1000 0.2513 0.001280
10000 0.2497 -0.000267
100000 0.2501 0.000107
 
p: 0.7000 theory: 0.2100 t: 100
n sim sim-theory
100 0.2171 0.007100
1000 0.2095 -0.000490
10000 0.2099 -0.000137
100000 0.2103 0.000321
 
p: 0.9000 theory: 0.0900 t: 100
n sim sim-theory
100 0.0999 0.009900
1000 0.0898 -0.000240
10000 0.0906 0.000568
100000 0.0908 0.000775</pre>
 
=={{header|Raku}}==
(formerly Perl 6)
<syntaxhighlight lang="raku" line>sub R($n, $p) { [+] ((rand < $p) xx $n).squish }
 
say 't= ', constant t = 100;
 
for .1, .3 ... .9 -> $p {
say "p= $p, K(p)= {$p*(1-$p)}";
for 10, 100, 1000 -> $n {
printf " R(%6d, p)= %f\n", $n, t R/ [+] R($n, $p)/$n xx t
}
}</syntaxhighlight>
{{out}}
<pre>t= 100
p= 0.1, K(p)= 0.09
R( 10, p)= 0.088000
R( 100, p)= 0.085600
R( 1000, p)= 0.089150
p= 0.3, K(p)= 0.21
R( 10, p)= 0.211000
R( 100, p)= 0.214600
R( 1000, p)= 0.211160
p= 0.5, K(p)= 0.25
R( 10, p)= 0.279000
R( 100, p)= 0.249200
R( 1000, p)= 0.250870
p= 0.7, K(p)= 0.21
R( 10, p)= 0.258000
R( 100, p)= 0.215400
R( 1000, p)= 0.209560
p= 0.9, K(p)= 0.09
R( 10, p)= 0.181000
R( 100, p)= 0.094500
R( 1000, p)= 0.091330</pre>
 
=={{header|Sidef}}==
{{trans|Perl 6Raku}}
<langsyntaxhighlight lang="ruby">func R(n,p) {
n.of { 1.rand < p ? 1 : 0}.sum;
}
Line 842 ⟶ 1,586:
printf (" R(n, p)= %f\n", t.of { R(n, p) }.sum/n / t);
}
}</langsyntaxhighlight>
 
{{out}}
Line 870 ⟶ 1,614:
 
=={{header|Tcl}}==
<langsyntaxhighlight lang="tcl">proc randomString {length probability} {
for {set s ""} {[string length $s] < $length} {} {
append s [expr {rand() < $probability}]
Line 902 ⟶ 1,646:
}
puts ""
}</langsyntaxhighlight>
{{out}}
<pre>
Line 924 ⟶ 1,668:
t=500, p=0.90, n= 4096, p(1-p)=0.090, sim=0.090, delta=0.08%
t=500, p=0.90, n=16384, p(1-p)=0.090, sim=0.090, delta=0.09%
</pre>
 
=={{header|Wren}}==
{{trans|Kotlin}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="wren">import "random" for Random
import "./fmt" for Fmt
 
var rand = Random.new()
var RAND_MAX = 32767
 
// just generate 0s and 1s without storing them
var runTest = Fn.new { |p, len, runs|
var cnt = 0
var thresh = (p * RAND_MAX).truncate
for (r in 0...runs) {
var x = 0
var i = len
while (i > 0) {
i = i - 1
var y = (rand.int(RAND_MAX + 1) < thresh) ? 1 : 0
if (x < y) cnt = cnt + 1
x = y
}
}
return cnt / runs / len
}
 
System.print("Running 1000 tests each:")
System.print(" p\t n\tK\tp(1-p)\t diff")
System.print("------------------------------------------------")
var fmt = "$.1f\t$6d\t$.4f\t$.4f\t$+.4f ($+.2f\%)"
for (ip in [1, 3, 5, 7, 9]) {
var p = ip / 10
var p1p = p * (1 - p)
var n = 100
while (n <= 1e5) {
var k = runTest.call(p, n, 1000)
Fmt.lprint(fmt, [p, n, k, p1p, k - p1p, (k - p1p) /p1p * 100])
n = n * 10
}
System.print()
}</syntaxhighlight>
 
{{out}}
Sample run:
<pre>
Running 1000 tests each:
p n K p(1-p) diff
------------------------------------------------
0.1 100 0.0919 0.0900 +0.0019 (+2.09%)
0.1 1000 0.0902 0.0900 +0.0002 (+0.23%)
0.1 10000 0.0900 0.0900 +0.0000 (+0.00%)
0.1 100000 0.0900 0.0900 -0.0000 (-0.03%)
 
0.3 100 0.2109 0.2100 +0.0009 (+0.43%)
0.3 1000 0.2102 0.2100 +0.0002 (+0.09%)
0.3 10000 0.2100 0.2100 -0.0000 (-0.01%)
0.3 100000 0.2100 0.2100 +0.0000 (+0.00%)
 
0.5 100 0.2527 0.2500 +0.0027 (+1.08%)
0.5 1000 0.2497 0.2500 -0.0003 (-0.10%)
0.5 10000 0.2501 0.2500 +0.0001 (+0.05%)
0.5 100000 0.2500 0.2500 -0.0000 (-0.01%)
 
0.7 100 0.2166 0.2100 +0.0066 (+3.14%)
0.7 1000 0.2107 0.2100 +0.0007 (+0.35%)
0.7 10000 0.2101 0.2100 +0.0001 (+0.05%)
0.7 100000 0.2100 0.2100 +0.0000 (+0.01%)
 
0.9 100 0.0970 0.0900 +0.0070 (+7.79%)
0.9 1000 0.0910 0.0900 +0.0010 (+1.14%)
0.9 10000 0.0901 0.0900 +0.0001 (+0.11%)
0.9 100000 0.0900 0.0900 +0.0000 (+0.04%)
</pre>
 
=={{header|zkl}}==
{{trans|C}}
<langsyntaxhighlight lang="zkl">fcn run_test(p,len,runs){
cnt:=0; do(runs){
pv:=0; do(len){
Line 937 ⟶ 1,755:
}
return(cnt.toFloat() / runs / len);
}</langsyntaxhighlight>
<langsyntaxhighlight lang="zkl">println("Running 1000 tests each:\n"
" p\t n\tK\tp(1-p)\t diff\n"
"-----------------------------------------------");
Line 950 ⟶ 1,768:
}
println();
}</langsyntaxhighlight>
{{out}}
<pre>
PureFox
9,476

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/234279...352185"