Percolation/Mean run density

Percolation/Mean run density is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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

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

The following vector has R10 = 3

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

Percolation theory states that

   K(p) = Rn / n = p(1 - p) as n tends to infinity.

Task

Any calculation of Rn / n for finite n is subject to randomnes so should be computed as the average of t runs, t >= 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

• s-Run on Wolfram mathworld.
• Percolation/Bond percolation
• Percolation/Site percolation

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%