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$, (and 0 is therefore $1-p$). Define a run of 1's as being a group of consecutive 1's in the vector bounded either by the limits of the vector or by a 0. Let the number of runs in a vector of length $n$ be $R_{n}$.

Percolation/Mean run density
You are encouraged to solve this task according to the task description, using any language you may know.

Percolation Simulation
This is a simulation of aspects of mathematical percolation theory.

For other percolation simulations, see Category:Percolation Simulations, or:
1D finite grid simulation
Mean run density
2D finite grid simulations
Site percolation | Bond percolation | Mean cluster density

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.

D

Translation of: python

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

      std.math;

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

auto newVect(in size_t len, in double prob) {

   return len.iota.map!(_ => uniform(0.0, 1.0) < prob).array;

}

size_t nRuns(R)(R vec) if (isForwardRange!R) {

   return vec.group.filter!q{ a[0] }.walkLength;

}

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

   return nRuns(newVect(n, prob)) / cast(double)n;

}

void main() {

   foreach (immutable p10; iota(1, 10, 2)) {
       immutable p = p10 / 10.0;
       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;
                           .reduce!q{a + b} / t;
           writefln("t=%3d, p=%4.2f, n=%5d, p(1-p)=%5.3f," ~
                    " sim=%5.3f, 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.090, sim=0.090, delta=0.5%
t=500, p=0.10, n= 4096, p(1-p)=0.090, sim=0.090, delta=0.2%
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.2%
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.250, delta=0.1%
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.210, delta=0.2%
t=500, p=0.70, n= 4096, p(1-p)=0.210, sim=0.210, delta=0.0%
t=500, p=0.70, n=16384, p(1-p)=0.210, sim=0.210, delta=0.1%

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

Tcl

<lang tcl>proc randomString {length probability} {

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

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

   }
   return $s

}

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

proc runs {str} {

   regexp -all {1+} $str

}

1. 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}]

}

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