Percolation/Mean run density

Let $v$ be a vector of $n$ values of either 1 or 0 where the probability of any value being 1 is $p$ ; the probability of a value being 0 is therefore $1-p$ . Define a run of 1s as being a group of consecutive 1s in the vector bounded either by the limits of the vector or by a 0. Let the number of such runs in a given vector of length $n$ be $R_{n}$ .

For example, the following vector has $R_{10}=3$

[1 1 0 0 0 1 0 1 1 1]
 ^^^       ^   ^^^^^

Percolation theory states that

K(p)=\lim _{n\to \infty }R_{n}/n=p(1-p)

Task

Any calculation of $R_{n}/n$ for finite $n$ is subject to randomness so should be computed as the average of $t$ runs, where $t\geq 100$ .

For values of $p$ of 0.1, 0.3, 0.5, 0.7, and 0.9, show the effect of varying $n$ on the accuracy of simulated $K(p)$ .

Show your output here.

See also

s-Run on Wolfram mathworld.

C

<lang c>#include <stdio.h>

include <stdlib.h>

// just generate 0s and 1s without storing them double run_test(double p, int len, int runs) { int r, x, y, i, cnt = 0, thresh = p * RAND_MAX;

for (r = 0; r < runs; r++) for (x = 0, i = len; i--; x = y) cnt += x < (y = rand() < thresh);

return (double)cnt / runs / len; }

int main(void) { double p, p1p, K; int ip, n;

puts( "running 1000 tests each:\n" " p\t n\tK\tp(1-p)\t diff\n" "-----------------------------------------------"); for (ip = 1; ip < 10; ip += 2) { p = ip / 10., p1p = p * (1 - p);

for (n = 100; n <= 100000; n *= 10) { K = run_test(p, n, 1000); printf("%.1f\t%6d\t%.4f\t%.4f\t%+.4f (%+.2f%%)\n", p, n, K, p1p, K - p1p, (K - p1p) / p1p * 100); } putchar('\n'); }

return 0; }</lang>

Output:

running 1000 tests each:
 p         n    K       p(1-p)       diff
-----------------------------------------------
0.1        100  0.0900  0.0900  -0.0001 (-0.06%)
0.1       1000  0.0899  0.0900  -0.0001 (-0.11%)
0.1      10000  0.0902  0.0900  +0.0002 (+0.17%)
0.1     100000  0.0900  0.0900  -0.0000 (-0.03%)

0.3        100  0.2110  0.2100  +0.0010 (+0.46%)
0.3       1000  0.2104  0.2100  +0.0004 (+0.19%)
0.3      10000  0.2100  0.2100  -0.0000 (-0.02%)
0.3     100000  0.2100  0.2100  -0.0000 (-0.01%)

0.5        100  0.2516  0.2500  +0.0016 (+0.66%)
0.5       1000  0.2498  0.2500  -0.0002 (-0.10%)
0.5      10000  0.2500  0.2500  +0.0000 (+0.01%)
0.5     100000  0.2500  0.2500  +0.0000 (+0.01%)

0.7        100  0.2162  0.2100  +0.0062 (+2.93%)
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.07%)
0.9       1000  0.0905  0.0900  +0.0005 (+0.57%)
0.9      10000  0.0901  0.0900  +0.0001 (+0.09%)
0.9     100000  0.0900  0.0900  +0.0000 (+0.03%)

C++

<lang cpp>#include <algorithm>

include <random>
include <vector>
include <iostream>
include <numeric>
include <iomanip>

using VecIt = std::vector<int>::const_iterator ;

//creates vector of length n, based on probability p for 1 std::vector<int> createVector( int n, double p ) {

  std::vector<int> result( n ) ;
  std::random_device rd ;
  std::mt19937 gen( rd( ) ) ;
  std::uniform_real_distribution<> dis( 0 , 1 ) ;
  for ( int i = 0 ; i < n ; i++ ) {
     double number = dis( gen ) ;
     if ( number <= p )

result[ i ] = 1 ;

     else

result[ i ] = 0 ;

  }
  return result ;

}

//find number of 1 runs in the vector int find_Runs( const std::vector<int> & numberVector ) {

  int runs = 0 ;
  VecIt found = numberVector.begin( ) ;
  while ( ( found = std::find( found , numberVector.end( ) , 1 ) )

!= numberVector.end( ) ) {

     runs++ ;
     while ( found != numberVector.end( ) && ( *found == 1 ) )

std::advance( found , 1 ) ;

     if ( found == numberVector.end( ) )

break ;

  }
  return runs ;

}

int main( ) {

  std::cout << "t = 100\n" ;
  std::vector<double> p_values { 0.1 , 0.3 , 0.5 , 0.7 , 0.9 } ;
  for ( double p : p_values ) { 
     std::cout << "p = " << p << " , K(p) = " << p * ( 1 - p ) << std::endl ;
     for ( int n = 10 ; n < 100000 ; n *= 10 ) {

std::vector<double> runsFound ; for ( int i = 0 ; i < 100 ; i++ ) { std::vector<int> ones_and_zeroes = createVector( n , p ) ; runsFound.push_back( find_Runs( ones_and_zeroes ) / static_cast<double>( n ) ) ; } double average = std::accumulate( runsFound.begin( ) , runsFound.end( ) , 0.0 ) / runsFound.size( ) ; std::cout << " R(" << std::setw( 6 ) << std::right << n << ", p) = " << average << std::endl ;

     }
  }
  return 0 ;

}</lang>

Output:

t = 100
p = 0.1 , K(p) = 0.09
  R(    10, p) = 0.088
  R(   100, p) = 0.0931
  R(  1000, p) = 0.09013
  R( 10000, p) = 0.089947
p = 0.3 , K(p) = 0.21
  R(    10, p) = 0.225
  R(   100, p) = 0.2089
  R(  1000, p) = 0.21043
  R( 10000, p) = 0.20991
p = 0.5 , K(p) = 0.25
  R(    10, p) = 0.271
  R(   100, p) = 0.253
  R(  1000, p) = 0.25039
  R( 10000, p) = 0.250278
p = 0.7 , K(p) = 0.21
  R(    10, p) = 0.264
  R(   100, p) = 0.2155
  R(  1000, p) = 0.20829
  R( 10000, p) = 0.209977
p = 0.9 , K(p) = 0.09
  R(    10, p) = 0.167
  R(   100, p) = 0.0928
  R(  1000, p) = 0.09071
  R( 10000, p) = 0.090341

D

Translation of: python

<lang d>import std.stdio, std.range, std.algorithm, std.random, std.math;

enum n = 100, p = 0.5, t = 500;

double meanRunDensity(in size_t n, in double prob) {

   return n.iota.map!(_ => uniform01 < prob)
          .array.uniq.sum / double(n);

}

void main() {

   foreach (immutable p; iota(0.1, 1.0, 0.2)) {
       immutable limit = p * (1 - p);
       writeln;
       foreach (immutable n2; iota(10, 16, 2)) {
           immutable n = 2 ^^ n2;
           immutable sim = t.iota.map!(_ => meanRunDensity(n, p))
                           .sum / t;
           writefln("t=%3d, p=%4.2f, n=%5d, p(1-p)=%5.5f, " ~
                    "sim=%5.5f, delta=%3.1f%%", t, p, n, limit, sim,
                    limit ? abs(sim - limit) / limit * 100 : sim*100);
       }
   }

}</lang>

Output:

t=500, p=0.10, n= 1024, p(1-p)=0.09000, sim=0.08949, delta=0.6%
t=500, p=0.10, n= 4096, p(1-p)=0.09000, sim=0.08976, delta=0.3%
t=500, p=0.10, n=16384, p(1-p)=0.09000, sim=0.08988, delta=0.1%

t=500, p=0.30, n= 1024, p(1-p)=0.21000, sim=0.20979, delta=0.1%
t=500, p=0.30, n= 4096, p(1-p)=0.21000, sim=0.21020, delta=0.1%
t=500, p=0.30, n=16384, p(1-p)=0.21000, sim=0.21005, delta=0.0%

t=500, p=0.50, n= 1024, p(1-p)=0.25000, sim=0.25016, delta=0.1%
t=500, p=0.50, n= 4096, p(1-p)=0.25000, sim=0.25026, delta=0.1%
t=500, p=0.50, n=16384, p(1-p)=0.25000, sim=0.24990, delta=0.0%

t=500, p=0.70, n= 1024, p(1-p)=0.21000, sim=0.21050, delta=0.2%
t=500, p=0.70, n= 4096, p(1-p)=0.21000, sim=0.20993, delta=0.0%
t=500, p=0.70, n=16384, p(1-p)=0.21000, sim=0.21009, delta=0.0%

t=500, p=0.90, n= 1024, p(1-p)=0.09000, sim=0.09019, delta=0.2%
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%

Fortran

<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

program percolation_mean_run_density

 implicit none
 integer :: i, p10, n2, n, t
 real :: p, limit, sim, delta
 data n,p,t/100,0.5,500/
 write(6,'(a3,a5,4a7)')'t','p','n','p(1-p)','sim','delta%'
 do p10=1,10,2
    p = p10/10.0
    limit = p*(1-p)
    write(6,'()')
    do n2=10,15,2
       n = 2**n2
       sim = 0
       do i=1,t
          sim = sim + mean_run_density(n,p)
       end do
       sim = sim/t
       if (limit /= 0) then
          delta = abs(sim-limit)/limit
       else
          delta = sim
       end if
       delta = delta * 100
       write(6,'(i3,f5.2,i7,2f7.3,f5.1)')t,p,n,limit,sim,delta
    end do
 end do

contains

 integer function runs(n, p)
   integer, intent(in) :: n
   real, intent(in) :: p
   real :: harvest
   logical :: q
   integer :: count, i
   count = 0
   q = .false.
   do i=1,n
      call random_number(harvest)
      if (harvest < p) then
         q = .true.
      else
         if (q) count = count+1
         q = .false.
      end if
   end do
   runs = count
 end function runs

 real function mean_run_density(n, p)
   integer, intent(in) :: n
   real, intent(in) :: p
   mean_run_density = real(runs(n,p))/real(n)
 end function mean_run_density

end program percolation_mean_run_density </lang>

$ ./f
  t    p      n p(1-p)    sim  delta%

500 0.10   1024  0.090  0.090  0.2
500 0.10   4096  0.090  0.090  0.2
500 0.10  16384  0.090  0.090  0.0

500 0.30   1024  0.210  0.210  0.2
500 0.30   4096  0.210  0.210  0.0
500 0.30  16384  0.210  0.210  0.0

500 0.50   1024  0.250  0.250  0.1
500 0.50   4096  0.250  0.250  0.1
500 0.50  16384  0.250  0.250  0.1

500 0.70   1024  0.210  0.210  0.1
500 0.70   4096  0.210  0.210  0.1
500 0.70  16384  0.210  0.210  0.0

500 0.90   1024  0.090  0.090  0.1
500 0.90   4096  0.090  0.090  0.4
500 0.90  16384  0.090  0.090  0.1

Haskell

<lang Haskell>import Control.Monad.Random import Control.Applicative import Text.Printf import Control.Monad import Data.Bits

data OneRun = OutRun | InRun deriving (Eq, Show)

randomList :: Int -> Double -> Rand StdGen [Int] randomList n p = take n . map f <$> getRandomRs (0,1)

 where f n = if (n > p) then 0 else 1

countRuns xs = fromIntegral . sum $

              zipWith (\x y -> x .&. xor y 1) xs (tail xs ++ [0])

calcK :: Int -> Double -> Rand StdGen Double calcK n p = (/ fromIntegral n) . countRuns <$> randomList n p

printKs :: StdGen -> Double -> IO () printKs g p = do

 printf "p= %.1f, K(p)= %.3f\n" p (p * (1 - p))
 forM_ [1..5] $ \n -> do
   let est = evalRand (calcK (10^n) p) g
   printf "n=%7d, estimated K(p)= %5.3f\n" (10^n::Int) est

main = do

 x <- newStdGen
 forM_ [0.1,0.3,0.5,0.7,0.9] $ printKs x

</lang>

./percolation
p= 0.1, K(p)= 0.090
n=     10, estimated K(p)= 0.000
n=    100, estimated K(p)= 0.130
n=   1000, estimated K(p)= 0.099
n=  10000, estimated K(p)= 0.090
n= 100000, estimated K(p)= 0.091
p= 0.3, K(p)= 0.210
n=     10, estimated K(p)= 0.200
n=    100, estimated K(p)= 0.250
n=   1000, estimated K(p)= 0.209
n=  10000, estimated K(p)= 0.209
n= 100000, estimated K(p)= 0.211
p= 0.5, K(p)= 0.250
n=     10, estimated K(p)= 0.200
n=    100, estimated K(p)= 0.290
n=   1000, estimated K(p)= 0.252
n=  10000, estimated K(p)= 0.250
n= 100000, estimated K(p)= 0.250
p= 0.7, K(p)= 0.210
n=     10, estimated K(p)= 0.300
n=    100, estimated K(p)= 0.200
n=   1000, estimated K(p)= 0.210
n=  10000, estimated K(p)= 0.209
n= 100000, estimated K(p)= 0.210
p= 0.9, K(p)= 0.090
n=     10, estimated K(p)= 0.200
n=    100, estimated K(p)= 0.090
n=   1000, estimated K(p)= 0.089
n=  10000, estimated K(p)= 0.095
n= 100000, estimated K(p)= 0.090

Icon and Unicon

The following works in both languages:

<lang unicon>procedure main(A)

   t := integer(A[2]) | 500

   write(left("p",8)," ",left("n",8)," ",left("p(1-p)",10)," ",left("SimK(p)",10))
   every (p := 0.1 | 0.3 | 0.5 | 0.7 | 0.9, n := 1000 | 2000 | 3000) do {
       Ka := 0.0
       every !t do {
           every (v := "", !n) do v ||:= |((?0.1 > p,"0")|"1")
           R := 0
           v ? while tab(upto('1')) do R +:= (tab(many('1')), 1)
           Ka +:= real(R)/n
           }
       write(left(p,8)," ",left(n,8)," ",left(p*(1-p),10)," ",left(Ka/t, 10))
       }

end</lang>

Output:

->pmrd
p        n        p(1-p)     SimK(p)   
0.1      1000     0.09000000 0.09021400
0.1      2000     0.09000000 0.08984799
0.1      3000     0.09000000 0.08993666
0.3      1000     0.21       0.21080999
0.3      2000     0.21       0.209953  
0.3      3000     0.21       0.210564  
0.5      1000     0.25       0.250024  
0.5      2000     0.25       0.25007399
0.5      3000     0.25       0.24975266
0.7      1000     0.21       0.21098799
0.7      2000     0.21       0.20987700
0.7      3000     0.21       0.21047333
0.9      1000     0.08999999 0.09016400
0.9      2000     0.08999999 0.09004800
0.9      3000     0.08999999 0.09023200
->

J

<lang J> NB. translation of python

NB. 'N P T' =: 100 0.5 500 NB. silliness

newv =: (> ?@(#&0))~ NB. generate a random binary vector. Use: N newv P runs =: {: + [: +/ 1 0&E. NB. add the tail to the sum of 1 0 occurrences Use: runs V mean_run_density =: [ %~ [: runs newv NB. perform experiment. Use: N mean_run_density P

main =: 3 : 0 NB.Usage: main T

T =. y
smoutput'  T  P    N    P(1-P) SIM   DELTA%'
for_P. 10 %~ >: +: i. 5 do.
  LIMIT =. (* -.) P
  smoutput 
  for_N. 2 ^ 10 + +: i. 3 do.
    SIM =. T %~ +/ (N mean_run_density P"_)^:(<T) 0
    smoutput 4 5j2 6 6j3 6j3 4j1 ": T, P, N, LIMIT, SIM, SIM (100 * [`(|@:(- % ]))@.(0 ~: ])) LIMIT
  end.
end.
EMPTY

) </lang> Session:

  main 500
  T  P    N    P(1-P) SIM   DELTA%

 500 0.10  1024 0.090 0.090 0.1
 500 0.10  4096 0.090 0.090 0.2
 500 0.10 16384 0.090 0.090 0.2

 500 0.30  1024 0.210 0.210 0.2
 500 0.30  4096 0.210 0.209 0.3
 500 0.30 16384 0.210 0.210 0.1

 500 0.50  1024 0.250 0.250 0.2
 500 0.50  4096 0.250 0.250 0.1
 500 0.50 16384 0.250 0.250 0.2

 500 0.70  1024 0.210 0.210 0.0
 500 0.70  4096 0.210 0.210 0.2
 500 0.70 16384 0.210 0.210 0.2

 500 0.90  1024 0.090 0.091 1.1
 500 0.90  4096 0.090 0.090 0.1
 500 0.90 16384 0.090 0.090 0.1

Perl

Translation of: Perl 6

<lang perl>sub R {

   my ($n, $p) = @_;
   my $r = join ,
   map { rand() < $p ? 1 : 0 } 1 .. $n;
   0+ $r =~ s/1+//g;

}

use constant t => 100;

printf "t= %d\n", t; for my $p (qw(.1 .3 .5 .7 .9)) {

   printf "p= %f, K(p)= %f\n", $p, $p*(1-$p);  
   for my $n (qw(10 100 1000)) {
       my $r; $r += R($n, $p) for 1 .. t; $r /= $n;
       printf " R(n, p)= %f\n", $r / t;
   }

}</lang>

Output:

t= 100
p= 0.100000, K(p)= 0.090000
 R(n, p)= 0.095000
 R(n, p)= 0.088100
 R(n, p)= 0.089420
p= 0.300000, K(p)= 0.210000
 R(n, p)= 0.225000
 R(n, p)= 0.208800
 R(n, p)= 0.210020
p= 0.500000, K(p)= 0.250000
 R(n, p)= 0.289000
 R(n, p)= 0.249900
 R(n, p)= 0.248980
p= 0.700000, K(p)= 0.210000
 R(n, p)= 0.262000
 R(n, p)= 0.213200
 R(n, p)= 0.209690
p= 0.900000, K(p)= 0.090000
 R(n, p)= 0.177000
 R(n, p)= 0.096200
 R(n, p)= 0.091730

Perl 6

<lang perl6>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

}</lang>

Output:

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

Python

<lang python>from __future__ import division from random import random from math import fsum

n, p, t = 100, 0.5, 500

def newv(n, p):

   return [int(random() < p) for i in range(n)]

def runs(v):

   return sum((a & ~b) for a, b in zip(v, v[1:] + [0]))

def mean_run_density(n, p):

   return runs(newv(n, p)) / n

for p10 in range(1, 10, 2):

   p = p10 / 10
   limit = p * (1 - p)
   print()
   for n2 in range(10, 16, 2):
       n = 2**n2
       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))</lang>

Output:

t=500 p=0.10 n= 1024 p(1-p)=0.090 sim=0.090 delta=0.2%
t=500 p=0.10 n= 4096 p(1-p)=0.090 sim=0.090 delta=0.0%
t=500 p=0.10 n=16384 p(1-p)=0.090 sim=0.090 delta=0.1%

t=500 p=0.30 n= 1024 p(1-p)=0.210 sim=0.210 delta=0.0%
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.3%
t=500 p=0.50 n= 4096 p(1-p)=0.250 sim=0.250 delta=0.0%
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.210 delta=0.0%
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=0.6%
t=500 p=0.90 n= 4096 p(1-p)=0.090 sim=0.090 delta=0.2%
t=500 p=0.90 n=16384 p(1-p)=0.090 sim=0.090 delta=0.0%

Racket

<lang racket>#lang racket (require racket/fixnum) (define t (make-parameter 100))

(define (Rn v)

 (define (inner-Rn rv idx b-1)
   (define b (fxvector-ref v idx))
   (define rv+ (if (and (= b 1) (= b-1 0)) (add1 rv) rv))
   (if (zero? idx) rv+ (inner-Rn rv+ (sub1 idx) b)))
 (inner-Rn 0 (sub1 (fxvector-length v)) 0))

(define ((make-random-bit-vector p) n)

 (for/fxvector
  #:length n ((i n))
  (if (<= (random) p) 1 0)))

(define (Rn/n l->p n) (/ (Rn (l->p n)) n))

(for ((p (in-list '(1/10 3/10 1/2 7/10 9/10))))

 (define l->p (make-random-bit-vector p))
 (define Kp (* p (- 1 p)))
 (printf "p = ~a\tK(p) =\t~a\t~a~%" p Kp (real->decimal-string Kp 4))
 (for ((n (in-list '(10 100 1000 10000))))
   (define sum-Rn/n (for/sum ((i (in-range (t)))) (Rn/n l->p n)))
   (define sum-Rn/n/t (/ sum-Rn/n (t)))
   (printf "mean(R_~a/~a) =\t~a\t~a~%"
           n n sum-Rn/n/t (real->decimal-string sum-Rn/n/t 4)))
 (newline))

(module+ test

 (require rackunit)
 (check-eq? (Rn (fxvector 1 1 0 0 0 1 0 1 1 1)) 3))</lang>

Output:

p = 1/10	K(p) =	9/100	0.0900
mean(R_10/10) =	3/40	0.0750
mean(R_100/100) =	221/2500	0.0884
mean(R_1000/1000) =	4469/50000	0.0894
mean(R_10000/10000) =	90313/1000000	0.0903

p = 3/10	K(p) =	21/100	0.2100
mean(R_10/10) =	231/1000	0.2310
mean(R_100/100) =	1049/5000	0.2098
mean(R_1000/1000) =	131/625	0.2096
mean(R_10000/10000) =	209873/1000000	0.2099

p = 1/2	K(p) =	1/4	0.2500
mean(R_10/10) =	297/1000	0.2970
mean(R_100/100) =	1263/5000	0.2526
mean(R_1000/1000) =	24893/100000	0.2489
mean(R_10000/10000) =	124963/500000	0.2499

p = 7/10	K(p) =	21/100	0.2100
mean(R_10/10) =	131/500	0.2620
mean(R_100/100) =	2147/10000	0.2147
mean(R_1000/1000) =	1049/5000	0.2098
mean(R_10000/10000) =	210453/1000000	0.2105

p = 9/10	K(p) =	9/100	0.0900
mean(R_10/10) =	169/1000	0.1690
mean(R_100/100) =	119/1250	0.0952
mean(R_1000/1000) =	4503/50000	0.0901
mean(R_10000/10000) =	89939/1000000	0.0899

Sidef

Translation of: Perl 6

<lang ruby>func R(n,p) {

   n.of { 1.rand < p ? 1 : 0}.sum;

}

say ('t=', const t = 100);

range(.1, .9, .2).each { |p|

   printf("p= %f, K(p)= %f\n", p, p*(1-p));
   [10, 100, 1000].each { |n|
       printf (" R(n, p)= %f\n", t.of { R(n, p) }.sum/n / t);
   }

}</lang>

Output:

t=100
p= 0.100000, K(p)= 0.090000
 R(n, p)= 0.099000
 R(n, p)= 0.105000
 R(n, p)= 0.099810
p= 0.300000, K(p)= 0.210000
 R(n, p)= 0.301000
 R(n, p)= 0.289800
 R(n, p)= 0.300720
p= 0.500000, K(p)= 0.250000
 R(n, p)= 0.481000
 R(n, p)= 0.501800
 R(n, p)= 0.498260
p= 0.700000, K(p)= 0.210000
 R(n, p)= 0.695000
 R(n, p)= 0.698400
 R(n, p)= 0.701220
p= 0.900000, K(p)= 0.090000
 R(n, p)= 0.910000
 R(n, p)= 0.898500
 R(n, p)= 0.899080

Tcl

<lang tcl>proc randomString {length probability} {

   for {set s ""} {[string length $s] < $length} {} {

append s [expr {rand() < $probability}]

   }
   return $s

}

By default, [regexp -all] gives the number of times that the RE matches

proc runs {str} {

   regexp -all {1+} $str

}

Compute the mean run density

proc mrd {t p n} {

   for {set i 0;set total 0.0} {$i < $t} {incr i} {

set run [randomString $n $p] set total [expr {$total + double([runs $run])/$n}]

   }
   return [expr {$total / $t}]

}

Parameter sweep with nested [foreach]

set runs 500 foreach p {0.10 0.30 0.50 0.70 0.90} {

   foreach n {1024 4096 16384} {

set theory [expr {$p * (1 - $p)}] set sim [mrd $runs $p $n] set diffpc [expr {abs($theory-$sim)*100/$theory}] puts [format "t=%d, p=%.2f, n=%5d, p(1-p)=%.3f, sim=%.3f, delta=%.2f%%" \ $runs $p $n $theory $sim $diffpc]

   }
   puts ""

}</lang>

Output:

t=500, p=0.10, n= 1024, p(1-p)=0.090, sim=0.090, delta=0.07%
t=500, p=0.10, n= 4096, p(1-p)=0.090, sim=0.090, delta=0.06%
t=500, p=0.10, n=16384, p(1-p)=0.090, sim=0.090, delta=0.17%

t=500, p=0.30, n= 1024, p(1-p)=0.210, sim=0.210, delta=0.23%
t=500, p=0.30, n= 4096, p(1-p)=0.210, sim=0.210, delta=0.09%
t=500, p=0.30, n=16384, p(1-p)=0.210, sim=0.210, delta=0.01%

t=500, p=0.50, n= 1024, p(1-p)=0.250, sim=0.250, delta=0.10%
t=500, p=0.50, n= 4096, p(1-p)=0.250, sim=0.250, delta=0.07%
t=500, p=0.50, n=16384, p(1-p)=0.250, sim=0.250, delta=0.08%

t=500, p=0.70, n= 1024, p(1-p)=0.210, sim=0.211, delta=0.33%
t=500, p=0.70, n= 4096, p(1-p)=0.210, sim=0.210, delta=0.00%
t=500, p=0.70, n=16384, p(1-p)=0.210, sim=0.210, delta=0.01%

t=500, p=0.90, n= 1024, p(1-p)=0.090, sim=0.091, delta=1.61%
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%

zkl

Translation of: C

<lang zkl>fcn run_test(p,len,runs){

  cnt:=0; do(runs){
     pv:=0; do(len){
        v:=0 + ((0.0).random(1.0)<p);  // 0 or 1, value of V[n]
        cnt += (pv<v);  // if v is 1 & prev v was zero, inc cnt
        pv = v;
     }
  }
  return(cnt.toFloat() / runs / len);

}</lang> <lang zkl>println("Running 1000 tests each:\n" " p\t n\tK\tp(1-p)\t diff\n" "-----------------------------------------------"); foreach p in ([0.1..0.9,0.2]) {

  p1p:=p*(1.0 - p);
  n:=100; while(n <= 100000) {
     K:=run_test(p, n, 1000);
     "%.1f\t%6d\t%.4f\t%.4f\t%+.4f (%+.2f%%)".fmt(

p, n, K, p1p, K - p1p, (K - p1p) / p1p * 100).println();

     n *= 10;
  }
  println();

}</lang>

Output:

Running 1000 tests each:
 p	   n	K	p(1-p)	     diff
-----------------------------------------------
0.1	   100	0.0903	0.0900	+0.0003 (+0.36%)
0.1	  1000	0.0900	0.0900	-0.0000 (-0.01%)
0.1	 10000	0.0901	0.0900	+0.0001 (+0.16%)
0.1	100000	0.0900	0.0900	+0.0000 (+0.01%)

0.3	   100	0.2115	0.2100	+0.0015 (+0.73%)
0.3	  1000	0.2105	0.2100	+0.0005 (+0.23%)
0.3	 10000	0.2098	0.2100	-0.0002 (-0.07%)
0.3	100000	0.2100	0.2100	+0.0000 (+0.00%)

0.5	   100	0.2521	0.2500	+0.0021 (+0.86%)
0.5	  1000	0.2503	0.2500	+0.0003 (+0.13%)
0.5	 10000	0.2500	0.2500	-0.0000 (-0.01%)
0.5	100000	0.2500	0.2500	-0.0000 (-0.00%)

0.7	   100	0.2151	0.2100	+0.0051 (+2.41%)
0.7	  1000	0.2103	0.2100	+0.0003 (+0.16%)
0.7	 10000	0.2100	0.2100	+0.0000 (+0.00%)
0.7	100000	0.2100	0.2100	-0.0000 (-0.01%)

0.9	   100	0.0979	0.0900	+0.0079 (+8.74%)
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.00%)