
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


{\displaystyle v}

 be a vector of 



n


{\displaystyle n}

 values of either 1 or 0 where the probability of any
value being 1 is 



p


{\displaystyle p}

; the probability of a value being 0 is therefore 



1
−
p


{\displaystyle 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


{\displaystyle n}

 be 




R

n




{\displaystyle R_{n}}

.

For example, the following vector has 




R

10


=
3


{\displaystyle R_{10}=3}



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

Percolation theory states that





K
(
p
)
=

lim

n
→
∞



R

n



/

n
=
p
(
1
−
p
)


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


Task
Any calculation of 




R

n



/

n


{\displaystyle R_{n}/n}

 for finite 



n


{\displaystyle n}

 is subject to randomness so should be
computed as the average of 



t


{\displaystyle t}

 runs, where 



t
≥
100


{\displaystyle t\geq 100}

.

For values of 



p


{\displaystyle p}

 of 0.1, 0.3, 0.5, 0.7, and 0.9, show the effect of varying 



n


{\displaystyle n}


on the accuracy of simulated 



K
(
p
)


{\displaystyle K(p)}

.

Show your output here.

See also
 s-Run on Wolfram mathworld.
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