Probability That a Drunk Person Can Cross the Bridge [gnuplot]

YouTube

Simulation [gnuplot]

Simulation model of bridge and river

Source code (PLT file)

reset

# Parameter --------------------
w    = 0.8            # Stride length
level = 5             # Level of drunk
ang = pi/3*level      # Movable range

N    = 5e4            # Number of trials
suc    = 0.           # Success (float: 0. = 0.0)
fail = 0.             # Failure
STOP = 40             # Number of plotting repeatedly

xd    = 0.            # Bridge
xu    = 20.
yd    = 0.
yu    = 4.
Lx    = 2. 
Ly    = 4.

# Setting --------------------
set term gif animate delay 4 size 1280,720
filename = sprintf("Lev.%d.gif", level)
set output filename

set nokey
unset grid
set size ratio ((yu+Ly)-(yd-Ly)) / ((xu+2*Lx)-(xd-2*Lx))

set xr[xd-2*Lx:xu+2*Lx]
set yr[yd-Ly:yu+Ly]
set xl "x" font 'Times New Roman:Italic, 26'
set yl "y" font 'Times New Roman:Italic, 26'
set xtics 2
set ytics 2
set tics font 'Times New Roman,24'

# Functions --------------------
N(i) = sprintf("{/Times:Italic N}  = %d\n", i)
P(i, j) = sprintf("{/Times:Italic p}  = %.5f\n", i/j)
Sf(sum, fail) = sprintf("{/Times:Italic s_f} = %.1f\n", sum/fail)

# Plot --------------------
# Bridge
set arrow 1 nohead from 0, yd-Ly to 0, yd lc -1 lw 4 back
set arrow 2 nohead from 0, yu to 0, yu+Ly lc -1 lw 4 back
set arrow 3 nohead from xu, yd-Ly to xu, yd lc -1 lw 4 back
set arrow 4 nohead from xu, yu to xu, yu+Ly lc -1 lw 4 back
set arrow 5 nohead from 0, 0 to xu, 0 lc -1 lw 4 back
set arrow 6 nohead from 0, yu to xu, yu lc -1 lw 4 back
set arrow 7 nohead from xd-Lx, yd-Ly to xd-Lx, yu+Ly lc 7 lw 1.5 back
set arrow 8 nohead from xu+Lx, yd-Ly to xu+Lx, yu+Ly lc 7 lw 1.5 back

# River
set object 1 rect from xd, yu to xu, yu+Ly back fc rgb "light-cyan" fs solid
set object 2 rect from xd, yd-Ly to xu, yd back fc rgb "light-cyan" fs solid

# Display level of drunk
set label 1 sprintf("Level %d", level) center at graph 0.5, 1.05 font 'TimesNewRoman, 28'

sum = 0.                      # Sum of sf (sf: steps on falling in a river)
dec = 0                       # Counter for decimation

do for [i = 1:N] {
    # Initiate value
    x1 = 0.
    y1 = yu/2.
    x2 = x1
    y2 = y1
    sf = 0                    # Steps on falling in a river
    num = 3                   # Number of object
        
    # Calculate digits of i
    buf = i
    digit = 0
    while(buf!=0){
        digit = digit + 1
        buf = buf / 10
    }

    # Start position
    set object num circle at x1, y1 size 0.08 fc rgb 'black' fs solid front

    # Display number of trials
    set label 2 left N(i) font 'Times New Roman, 28' at graph 0.45, 0.94
    
    # Draw initiate state for STOP frame(i=1)
    if(i==1){
        do for[j=1:STOP]{ plot 1/0 }
    }
    
    # Loop
    while( (x2>=xd && x2<xu && y2>yd && y2<yu)||(x2<xd && x2>xd-Lx)||(x2>=xu && x2<xu+Lx) ) {
        sf = sf + 1
        num = num + 1
        x1 = x2                         # Old position
        y1 = y2
        r = rand(0)
        x2 = x2 + w*cos(ang*(r-0.5))    # New position
        y2 = y2 + w*sin(ang*(r-0.5))
        
        # Display a point and a line
        set object num circle at x2, y2 size 0.08 fc rgb 'black' fs solid front
        set arrow num+7 nohead from x1, y1 to x2, y2 lc -1 lw 2

        # SUCCESS: Reach the goal
        if(x2>=xu+Lx){
            suc = suc + 1
            set label 3 P(suc, i) left font 'Times New Roman, 28' at graph 0.45, 0.87
            
            if(fail!=0){
                set label 4 Sf(sum, fail) left font 'Times New Roman, 28' at graph 0.45, 0.80
            }
            
            # Turn path red
            t = 3
            set object t fc rgb 'red'
            while(t<num){
                t = t+1
                set object t fc rgb 'red'
                set arrow t+7 lc 7
            }
        }
        
        # FAILURE: Fall in a river 
        if ( (x2>=xd && x2<=xu && y2<=yd) || (x2>=xd && x2<=xu && y2>=yu) || x2<=xd-Lx){
            fail = fail + 1
            sum = sum + sf    # Sum of sf (sf: steps on falling in a river)
            
            set label 3 P(suc, i) left font 'Times New Roman, 28' at graph 0.45, 0.87
            set label 4 Sf(sum, fail) left font 'Times New Roman, 28' at graph 0.45, 0.80
        }
        
        # Plot process in a trial
        # i = 1 -> 10
        if(i<=10){
            plot 1/0
        }
        # i = 11 -> 20
        if(i>=11 && i<=20){
            dec = dec + 1
            if( (dec%3==0) || (x2<xd) || (x2>=xu) || (y2<=yd) || (y2>=yu) ){
                plot 1/0        # Decimate by 3 steps (from Cond.1)
            }
        }
    }
    
    # Plot results of each trial (i = 20 -> N)
    if(i>20 && i>=10**(digit-1) && i<10**digit && i%(10**(digit-1)/5)==0){
        plot 1/0                # Decimate
    }
    
    # Draw final state for STOP steps(i=N)
    if(i==N){
        do for[j=1:STOP]{ plot 1/0 }
    }
    
    # Remove the path for preparing the next trial
    while(num>=4){
        unset object num
        unset arrow num+7
        num = num - 1    
    }
}

set out

Output (Lev.1.gif → Lev.6.gif)

Extra : Graph representing change in the success probability and $s_f$

Source code (PLT file)

reset
set term png size 1000, 600
set output "Graph.png"

set multiplot layout 1, 2 title "{/Times:Italic N}=50000" font 'Times New Roman, 20'
set size ratio 1.3
set nokey
set grid
set xr[1:6]
set xtics 1
set xl 'Level' font 'Times New Roman:Normal, 22'

# Left side
set yr[0:1]
set ytics 0.1
set yl 'Success probability' font 'Times New Roman:Normal, 22'
plot '-' w linesp lt 1 lw 2 lc -1 ps 1 pt 7
    1, 0.84492
    2, 0.30508
    3, 0.06662
    4, 0.00752
    5, 0.00012
    6, 0.00000
    e

# Right side
set yr[0:20]
set ytics 2
set yl "Average numbers of steps  {/Times:Italic=24 s_f}" font 'Times New Roman, 22
plot '-' w linesp lt 1 lw 2 lc 3 ps 1 pt 7
    1, 19.0
    2, 15.2
    3, 13.8
    4, 13.5
    5, 15.0
    6, 14.5
    e

unset multiplot
set out

Output