Return to JUNTO

JUNTO Practice: Probability, Connected Components

Discussed on August 25, 2020.

From "Twenty problems in probability":

An m by n checkerboard is colored randomly: each square is randomly painted white or black. We say that two squares, p and q, are in the same connected monochromatic component (or component, in short) if there is a sequence of squares, all of the same color, starting at p and ending at q, in which succesive squares in the sequence share a common side.

You have two goals:

Please provide an answer and a justification.

NOTE: The goal of this problem is different than the goal of the original problem from "Twenty problems in probability".

Solutions

Click to see:

Oscar Martinez

JUNTO_components.py

import numpy as np
from scipy.ndimage.measurements import label
import pandas as pd
from numba import jit


def generate_board(m=8, n=8):
    return np.random.randint(0, 2, (m, n))


def generate_boards(m=8, n=8, simulations=100000):
    return np.random.randint(0, 2, (simulations, m, n))


def count_components_multidimensional(boards):
    structure = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
    components = [
        label(b, structure)[1] + label(1 - b, structure)[1]
        for b in boards
    ]
    return components


def count_components(board):
    structure = np.array([[0, 1, 0], [1, 1, 1], [0, 1, 0]])
    _, ncomponents_black = label(board, structure)
    _, ncomponents_white = label(1 - board, structure)
    return ncomponents_white + ncomponents_black


def simulate(m=8, n=8, simulations=100000):
    return np.mean(
        [
            count_components(generate_board(m, n))
            for _ in range(simulations)
        ]
    )


def simulate_multidimensional(m=8, n=8, simulations=100000):
    return np.mean(
        count_components_multidimensional(
            generate_boards(m, n, simulations)
        )
    )


@jit
def generate_dataset(max_m, max_n, reps=10):
    m_list = []
    n_list = []
    components_list = []
    for m in range(1, max_m):
        for n in range(1, max_n):
            for i in range(reps):
                components = count_components(
                    generate_board(m, n)
                )
                m_list.append(m)
                n_list.append(n)
                components_list.append(components)

    return pd.DataFrame(
        {
            "m": m_list,
            "n": n_list,
            "components": components_list,
        }
    )


@jit
def generate_mean_dataset(max_m, max_n, reps=10):
    m_list = []
    n_list = []
    components_list = []
    for m in range(1, max_m):
        for n in range(1, max_n):
            components_samples = []
            for i in range(reps):
                components = count_components(
                    generate_board(m, n)
                )
                components_samples.append(components)
                components_mean = np.mean(
                    components_samples
                )
            m_list.append(m)
            n_list.append(n)
            components_list.append(components_mean)

    return pd.DataFrame(
        {
            "m": m_list,
            "n": n_list,
            "components": components_list,
        }
    )

# We import the functions from the JUNTO_components module that solved part (a) of the problem
from JUNTO_components import *
from mpl_toolkits import mplot3d
# Initial imports
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
# Generate/Save/Load dataset of mean components

#df = generate_mean_dataset(50,50,1000)
#df.to_csv("JUNTO_components_mean_dataset.csv")
df = pd.read_csv("JUNTO_components_mean_dataset.csv")
df.head(10)
Unnamed: 0 m n components mn
0 0 1 1 1.000 1.000000
1 1 1 2 1.487 0.500000
2 2 1 3 1.988 0.333333
3 3 1 4 2.476 0.250000
4 4 1 5 2.992 0.200000
5 5 1 6 3.490 0.166667
6 6 1 7 4.017 0.142857
7 7 1 8 4.449 0.125000
8 8 1 9 4.965 0.111111
9 9 1 10 5.509 0.100000
# Generate/Save/Load dataset of sample component counts

#df_samples = generate_dataset(100,100,100)
#df_samples.to_csv("JUNTO_components_samples_dataset.csv")
df_samples = pd.read_csv("JUNTO_components_samples_dataset.csv")
df_samples.head(10)
Unnamed: 0 m n components
0 0 1 1 1
1 1 1 1 1
2 2 1 1 1
3 3 1 1 1
4 4 1 1 1
5 5 1 1 1
6 6 1 1 1
7 7 1 1 1
8 8 1 1 1
9 9 1 1 1
# matplotlib of samples data
ax = plt.axes(projection='3d')
ax.plot3D(df_samples.m,df_samples.n,df_samples.components)
[<mpl_toolkits.mplot3d.art3d.Line3D at 0x11c799fd0>]
# matplotlib of mean data
ax = plt.axes(projection='3d')
ax.plot3D(df.m,df.n,df.components)
[<mpl_toolkits.mplot3d.art3d.Line3D at 0x11f11f190>]
# We import interactive plotting data for 3D plots
import plotly
import plotly.graph_objs as go
plotly.offline.init_notebook_mode()
def interactive_plot(x,y,z):
    '''
    Creates an interactive plot from 3 iterables of values
    '''
    trace = go.Scatter3d(
    x=x,
    y=y,
    z=z,
    mode="markers",
    marker={'size':10,'opacity':0.8}
    )

    layout = go.Layout(
        margin={'l': 0, 'r': 0, 'b': 0, 't': 0}
    )

    data = [trace]

    plot_figure = go.Figure(data = data, layout = layout)

    plotly.offline.iplot(plot_figure)
    
    return True    
# mean data interactive plot
interactive_plot(
    x=df.m,
    y=df.n,
    z=df.components)
True
# import sklearn and numpy for regression to estimate function
from sklearn import linear_model
import numpy as np
# we try fitting a simple linear regression
reg = linear_model.LinearRegression().fit(df_samples[["m","n"]],df_samples["components"])
print(reg.score(df[["m","n"]],df["components"]))
print(reg.coef_)
-0.9583397579347497
[6.93436765 6.93503388]
reg.predict(np.array([1,1]).reshape(1,-1))
array([-314.73500451])
# fit a Generalized Linear Model with Exponential distribution to better reflect the concavity of data
glm_reg = linear_model.GammaRegressor(alpha=1,max_iter=1000)
glm_reg.fit(df_samples[["m","n"]],df_samples["components"])
GammaRegressor(alpha=1, max_iter=1000)
print(glm_reg.coef_)
print(glm_reg.intercept_)
print(np.exp(glm_reg.intercept_ + 2*glm_reg.coef_[0] + 2*glm_reg.coef_[1]))
print(glm_reg.__dict__)
[0.0242958 0.0242983]
3.083342365164521
24.059515135973623
{'alpha': 1, 'fit_intercept': True, 'link': 'log', 'solver': 'lbfgs', 'max_iter': 1000, 'tol': 0.0001, 'warm_start': False, 'verbose': 0, '_family_instance': <sklearn._loss.glm_distribution.GammaDistribution object at 0x12805eb90>, '_link_instance': <sklearn.linear_model._glm.link.LogLink object at 0x12805ec50>, 'n_iter_': 23, 'intercept_': 3.083342365164521, 'coef_': array([0.0242958, 0.0242983])}

E(Ci) = exp(2.0396002689775847 + mi * 0.04474463 + ni * 0.0446625)

glm_reg.score(df[["m","n"]],df["components"])
0.6909683175752671
glm_reg.predict(np.array([8,8]).reshape(1,-1))
array([32.2041439])
# generate simulated data with exponential glm predictor to see how similar it looks to the truth
m_sims = []
n_sims = []
predictions = []           
for m in range(1,50):
    for n in range(1,50):
        prediction = glm_reg.predict(np.array([m,n]).reshape(1,-1))
        m_sims.append(m)
        n_sims.append(n)
        predictions.append(prediction[0])
interactive_plot(m_sims,n_sims,predictions)
True
# we try approximating the function with a neural network instead
import torch.nn as nn
import torch.nn.functional as F
import torch
# Version 1 of the neural network we will ultimately use
# 1 Hidden layer, 2 inputs, 1 output, tanh activation
class ComponentNet(nn.Module):
    def __init__(self):
        super().__init__()
        self.fc1 = nn.Linear(2,4)
        self.fc2 = nn.Linear(4,1)
        
    def forward(self, x):
        out = torch.tanh(self.fc1(x))
        out = self.fc2(out)
        
        return out
# we define a function to use our model (for demo purposes, use this instead of retraining)
class MeanComponentsInference():
    def __init__(self):
        self.trained_model = ComponentNet()
        self.trained_model.load_state_dict(torch.load("component_net.pt"))
        self.trained_model.eval()
        
    def infer(self,m,n):
        return self.trained_model(torch.tensor([m,n]).float())[0] - 11.7421
# torch acceptable format
X = torch.from_numpy(df_samples[["m","n"]].values)
y = torch.from_numpy(df_samples[["components"]].values)
# train and test split
n_samples = X.shape[0]
n_val = int(0.2 * n_samples)

shuffled_indices = torch.randperm(n_samples)

train_indices = shuffled_indices[:-n_val]
val_indices = shuffled_indices[-n_val:]

print(train_indices, val_indices)

train_X = X[train_indices]
val_X = X[val_indices]

train_y = y[train_indices]
val_y = y[val_indices]
tensor([589373,  45532, 282343,  ..., 640957, 148432, 450993]) tensor([736863,  91208, 142337,  ..., 321840,  10220, 171306])
# our data is large and will crash kernel unless it is batched with data loaders
train_loader = torch.utils.data.DataLoader(torch.utils.data.TensorDataset(train_X,train_y),batch_size=2000,shuffle=True)
val_loader = torch.utils.data.DataLoader(torch.utils.data.TensorDataset(val_X,val_y),batch_size=2000,shuffle=True)
for batch in train_loader:
    print(batch)
    break
[tensor([[86, 13],
        [51, 70],
        [50, 58],
        ...,
        [21, 87],
        [93, 82],
        [40, 64]]), tensor([[ 179],
        [ 552],
        [ 428],
        ...,
        [ 283],
        [1050],
        [ 407]])]
def training_loop(n_epochs, optimizer, model, loss_fn, train_X,train_y,test_X,test_y):
    # vanilla training loop for small datasets
    for epoch in range(1,n_epochs+1):
        output = model(train_X)
        loss_train = loss_fn(output,train_y.unsqueeze(1))
        
        val_output = model(test_X)
        loss_val = loss_fn(val_output,test_y.unsqueeze(1))
        
        optimizer.zero_grad()
        loss_train.backward()
        optimizer.step()
        
        if epoch == 1 or epoch % 5 ==0:
             print(f"Epoch {epoch}, Training loss {loss_train.item():.4f},"
            f" Validation loss {loss_val.item():.4f}")
def training_loop_loader(n_epochs, optimizer, model, loss_fn, train_loader,val_loader):
    # training loop that uses batches to spare the computer
    for epoch in range(1,n_epochs+1):
        for train_X,train_y in train_loader:
            batch_size = train_X.shape[0]
            outputs = model(train_X.float())
            loss_train = loss_fn(outputs,train_y.float())
            
            optimizer.zero_grad()
            loss_train.backward()
            optimizer.step()
        
        if epoch == 1 or epoch % 5 ==0:
             print(f"Epoch {epoch}, Training loss {loss_train.item():.4f}")
# we instantiate the model and use Adam optimizer for learning rate optimization, use mean squared error for our loss
import torch.optim as optim

model = ComponentNet()
optimizer = optim.Adam(model.parameters(),lr=1e-2)
loss_fn = nn.MSELoss()
# train the model
training_loop_loader(
    n_epochs=500,
    optimizer=optimizer,
    model = model,
    loss_fn = loss_fn,
    train_loader = train_loader,
    val_loader=val_loader,
)
Epoch 1, Training loss 1537.1517
Epoch 5, Training loss 875.4327
Epoch 10, Training loss 1122.5061
Epoch 15, Training loss 1317.2239
Epoch 20, Training loss 480.9785
Epoch 25, Training loss 809.4860
Epoch 30, Training loss 617.5715
Epoch 35, Training loss 674.5927
Epoch 40, Training loss 724.4601
Epoch 45, Training loss 932.1520
Epoch 50, Training loss 650.9850
Epoch 55, Training loss 480.2036
Epoch 60, Training loss 562.9608
Epoch 65, Training loss 2184.5410
Epoch 70, Training loss 1014.2282
Epoch 75, Training loss 969.5809
Epoch 80, Training loss 785.1955
Epoch 85, Training loss 735.7379
Epoch 90, Training loss 567.9627
Epoch 95, Training loss 473.5512
Epoch 100, Training loss 626.4579
Epoch 105, Training loss 812.8909
Epoch 110, Training loss 476.4218
Epoch 115, Training loss 1138.5979
Epoch 120, Training loss 395.0683
Epoch 125, Training loss 487.9115
Epoch 130, Training loss 1115.0310
Epoch 135, Training loss 633.8607
Epoch 140, Training loss 657.7079
Epoch 145, Training loss 965.3415
Epoch 150, Training loss 1497.4521
Epoch 155, Training loss 452.6031
Epoch 160, Training loss 773.6358
Epoch 165, Training loss 567.7449
Epoch 170, Training loss 645.7665
Epoch 175, Training loss 437.2297
Epoch 180, Training loss 415.5421
Epoch 185, Training loss 1214.3215
Epoch 190, Training loss 852.4963
Epoch 195, Training loss 844.5226
Epoch 200, Training loss 455.6103
Epoch 205, Training loss 768.0081
Epoch 210, Training loss 373.9290
Epoch 215, Training loss 357.3707
Epoch 220, Training loss 468.6095
Epoch 225, Training loss 623.4183
Epoch 230, Training loss 1438.2423
Epoch 235, Training loss 920.7835
Epoch 240, Training loss 713.5925
Epoch 245, Training loss 1868.9651
Epoch 250, Training loss 899.1036
Epoch 255, Training loss 547.0290
Epoch 260, Training loss 618.0707
Epoch 265, Training loss 1257.9821
Epoch 270, Training loss 487.6640
Epoch 275, Training loss 611.5565
Epoch 280, Training loss 528.6788
Epoch 285, Training loss 358.4393
Epoch 290, Training loss 867.3081
Epoch 295, Training loss 654.2693
Epoch 300, Training loss 499.9281
Epoch 305, Training loss 648.3791
Epoch 310, Training loss 639.1547
Epoch 315, Training loss 519.8954
Epoch 320, Training loss 752.7217
Epoch 325, Training loss 915.1215
Epoch 330, Training loss 576.0844
Epoch 335, Training loss 970.7849
Epoch 340, Training loss 543.1246
Epoch 345, Training loss 329.6135
Epoch 350, Training loss 577.3944
Epoch 355, Training loss 781.9626
Epoch 360, Training loss 699.2152
Epoch 365, Training loss 712.3376
Epoch 370, Training loss 865.6357
Epoch 375, Training loss 2515.6562
Epoch 380, Training loss 815.8028
Epoch 385, Training loss 531.0990
Epoch 390, Training loss 541.7260
Epoch 395, Training loss 726.5151
Epoch 400, Training loss 530.4302
# get the validation loss
with torch.no_grad():
    for val_X, val_y in val_loader:
        batch_size = val_X.shape[0]
        outputs = model(val_X.float())
        loss = loss_fn(outputs, val_y.float())
        print(f"Val loss {loss}")
Val loss 218246.59375
Val loss 223935.046875
Val loss 217332.921875
Val loss 223835.984375
Val loss 224207.796875
Val loss 229332.65625
Val loss 222901.25
Val loss 218963.515625
Val loss 219126.078125
Val loss 211241.25
Val loss 209158.34375
Val loss 222849.84375
Val loss 221020.390625
Val loss 217618.9375
Val loss 226837.078125
Val loss 230126.65625
Val loss 219834.03125
Val loss 213531.640625
Val loss 226420.8125
Val loss 217937.015625
Val loss 231883.734375
Val loss 222083.328125
Val loss 224181.859375
Val loss 221517.796875
Val loss 224796.953125
Val loss 220592.765625
Val loss 220406.453125
Val loss 219560.1875
Val loss 225369.265625
Val loss 232987.3125
Val loss 233223.59375
Val loss 214066.03125
Val loss 241962.453125
Val loss 217781.71875
Val loss 223730.28125
Val loss 215474.84375
Val loss 225510.96875
Val loss 218529.59375
Val loss 221289.828125
Val loss 215430.453125
Val loss 221082.484375
Val loss 238636.265625
Val loss 228500.953125
Val loss 232253.375
Val loss 218949.640625
Val loss 243541.890625
Val loss 226871.90625
Val loss 219967.109375
Val loss 216208.65625
Val loss 226055.15625
Val loss 228057.390625
Val loss 236820.890625
Val loss 222420.640625
Val loss 232826.703125
Val loss 232815.78125
Val loss 226873.328125
Val loss 218559.921875
Val loss 220140.59375
Val loss 216617.84375
Val loss 218102.875
Val loss 221078.84375
Val loss 221333.90625
Val loss 225427.375
Val loss 235122.609375
Val loss 211991.59375
Val loss 225973.390625
Val loss 217390.46875
Val loss 234201.71875
Val loss 219938.125
Val loss 218693.78125
Val loss 228797.46875
Val loss 231421.53125
Val loss 215438.375
Val loss 227202.3125
Val loss 229014.859375
Val loss 220406.75
Val loss 224926.1875
Val loss 238787.46875
Val loss 223827.8125
Val loss 236367.5625
Val loss 222431.359375
Val loss 230030.859375
Val loss 222160.453125
Val loss 239669.390625
Val loss 217799.046875
Val loss 221047.5625
Val loss 226343.65625
Val loss 218722.890625
Val loss 236069.296875
Val loss 235222.90625
Val loss 237637.125
Val loss 209434.65625
Val loss 225404.0
Val loss 228975.1875
Val loss 225586.140625
Val loss 215879.140625
Val loss 232287.53125
Val loss 228649.359375
Val loss 134665.484375
# save our model for later use
#torch.save(model.state_dict(),"component_net.pt")
# generate predictions to plot our model against the data
component_net = MeanComponentsInference()
with torch.no_grad():
    m_sims = []
    n_sims = []
    predictions = []
    for m in range(1,50):
        for n in range(1,50):
            prediction = component_net.infer(m,n)
            m_sims.append(m)
            n_sims.append(n)
            predictions.append(prediction)
# theres a weird bug in the model where we need to subtract by 12 to get better predictions, we also unpack the tensor
modded_preds = [p - 12 for p in predictions]
modded_preds = [p[0] for p in modded_preds]
---------------------------------------------------------------------------

IndexError                                Traceback (most recent call last)

<ipython-input-72-072f154faeb5> in <module>
      1 modded_preds = [p - 12 for p in predictions]
----> 2 modded_preds = [p[0] for p in modded_preds]


<ipython-input-72-072f154faeb5> in <listcomp>(.0)
      1 modded_preds = [p - 12 for p in predictions]
----> 2 modded_preds = [p[0] for p in modded_preds]


IndexError: invalid index to scalar variable.
interactive_plot(x=m_sims,
    y=n_sims,
    z=predictions)
True
# eyeball our predictions
for i,n in enumerate(n_sims):
    print(f"n:{n} m:{m_sims[i]} prediction:{predictions[i]}")
    if i>10:
        break
n:1 m:1 prediction:1.3522148132324219
n:2 m:1 prediction:2.075206756591797
n:3 m:1 prediction:2.769237518310547
n:4 m:1 prediction:3.435466766357422
n:5 m:1 prediction:4.074779510498047
n:6 m:1 prediction:4.688121795654297
n:7 m:1 prediction:5.276439666748047
n:8 m:1 prediction:5.840404510498047
n:9 m:1 prediction:6.380809783935547
n:10 m:1 prediction:6.898204803466797
n:11 m:1 prediction:7.393321990966797
n:12 m:1 prediction:7.866588592529297
# we see that our predictions are very closely in line with the means
component_net = MeanComponentsInference()
df["nn_pred"] = df.apply(lambda r: component_net.infer(r.m,r.n).item(),axis=1)
np.mean(np.square(df["components"] - df["nn_pred"]))
df["components"].plot()
df["nn_pred"].plot()
<matplotlib.axes._subplots.AxesSubplot at 0x128e03710>

John Lekberg

I believe that the expected number of connected components for an 8 by 8 checkerboard is

E(8, 8) ≈ 14.3

I believe that the general formula for the expected number of connected components in a n by m checkerboard is

E(n, m) ≈ 0.16 + 0.36 n + 0.36 m + 0.13 n m


To compute the expected number of connected components for an n by m checkerboard, I created a Python function, expected_connected_components:

import numpy as np
import scipy.ndimage as ndi


def expected_connected_components(*, dim, n_samples):
    """Compute the expected number of connected components
    on a checkerboard with given dimensions.
    
    dim -- Tuple[int, int]. Checkerboard dimensions. E.g. (8, 8)
    n_samples -- int. Number of samples to take. E.g. 10_000.
    """
    n, m = dim
    white_hot = np.random.randint(0, 2, (n_samples, n, m))
    black_hot = 1 - white_hot

    label_structure = np.array(
        [
            [[0, 0, 0], [0, 0, 0], [0, 0, 0]],
            [[0, 1, 0], [1, 1, 1], [0, 1, 0]],
            [[0, 0, 0], [0, 0, 0], [0, 0, 0]],
        ]
    )

    _, total_cc_white = ndi.label(
        white_hot, structure=label_structure
    )
    _, total_cc_black = ndi.label(
        black_hot, structure=label_structure
    )
    total_cc = total_cc_white + total_cc_black
    expected_cc = total_cc / n_samples

    return expected_cc


expected_connected_components(dim=(8, 8), n_samples=100_000)
14.33688

To develop a general formula in terms of n and m, I created a Python function, generate_cc_dataset, that generates a DataFrame of expectation data:

import pandas as pd


def generate_cc_dataset(*, dim_domain, n_samples):
    """Generate a DataFrame of expected connected
    components for different dimensions.
    
    dim_domain -- Iterable[Tuple[int, int]]. Different
        dimensions to collect data for. E.g. [(8, 8), (3, 5)].
    n_samples -- int. Number of samples to take. See
        `expected_connected_components` for details. E.g. 10_000.
    
    Output: DataFrame["n", "m", "cc"]
    """
    records = []
    for dim in dim_domain:
        n, m = dim
        cc = expected_connected_components(
            dim=dim, n_samples=n_samples
        )
        records.extend(
            [
                {"n": n, "m": m, "cc": cc},
                {"n": m, "m": n, "cc": cc},
            ]
        )
    return pd.DataFrame(records).drop_duplicates()

I generate a training dataset for all n and m such that:

from itertools import combinations_with_replacement

d_train = generate_cc_dataset(
    dim_domain=combinations_with_replacement(
        range(1, 33), 2
    ),
    n_samples=10000,
)
d_train.sample(n=5)
       n   m        cc
101   20   2   13.3429
1040  29  31  139.9753
19    10   1    5.4940
840   18  30   88.6639
387   23   7   32.0417

And I generate a testing dataset for a some n and m such that:

from random import randint

d_test = generate_cc_dataset(
    dim_domain=[
        (randint(100, 500), randint(100, 500))
        for _ in range(50)
    ],
    n_samples=1000,
)
d_test.sample(n=5)
      n    m         cc
78  464  141   8821.674
52  241  144   4701.318
9   379  343  17367.269
50  143  319   6172.145
63  276  365  13488.091

Now I need to develop a model for the expected number of connected components. For an n by m checkerboard,

As a result, I believe that I can model the expected number of connected components as a polynomial like this:

E ≈ C1 + C2n + C3m + C4nm

Here's a Python function, model_X_y that splits my training (and testing) data into

def model_X_y(dataset):
    """
    Split a dataset into X and y based on the model
    
    E ~ 1 + n + m + nm
    
    Input: DataFrame["n", "m", "cc"]
    
    Output:
    
    - X: DataFrame["n", "m", "nm"]
    - y: Series["cc"]
    """
    X = dataset[["n", "m"]]
    X["nm"] = X["n"] * X["m"]
    y = dataset[["cc"]]
    return X, y

I train my model on the training data:

from sklearn.linear_model import LinearRegression

model = LinearRegression()
train_X, train_y = model_X_y(d_train)
model.fit(train_X, train_y)

Then I test the predictive ability of the model on the test dataset:

from sklearn.metrics import r2_score

test_X, test_y = model_X_y(d_test)
test_y_pred = model.predict(test_X)

r2_score(test_y, test_y_pred)
0.9999898697943694

The predictive r2 score is 99.999%. As a result, I think this model does a good job predicting the estimated number of connected components. Here are the model intercept and coefficients:

model.intercept_
array([0.15607795])
model.coef_
array([[0.3608306 , 0.3608306 , 0.13134395]])

As a result, I believe that the formula for predicting the expected number of connected components is:

E(n, m) ≈ 0.16 + 0.36 n + 0.36 m + 0.13 n m

Daniel Bassett

Part 1: Expected(components) = ~11

I arrived at 11 expected components given an 8x8 board. I ran a simulation using boolean arrays to simulate a "black square or white square". While not exactly accurate, my answer does fall within the expected range of mn/8<E(C)<(m+2)(n+2)/6, given from the original Putnam question.

Part 2: Due to time constraints, I did not answer the prompt exactly, but rather proved the expected range of connected components given an MxN board. See the comments in the bottom of my solution for this proof.


import numpy as np

comp = []


def checkerboard():
    """Creates a 'checkerboard' in the form of an 8 by 8
    binary matrix. '1' signifies a black square, '0' a white
    square.
    """
    Matrix = np.random.randint(2, size=(8, 8))

    return Matrix


def test(board):
    """Simulates the connected components on the board. 
    Iterating through the board, squares are connected
    only if two true binaries are next to one another.
    That is, they are connected via column or row. I use 
    a negated conjunction to ensure there is not double 
    counting of components. This logic does not cover the
    breadth of possibilities, but for the sake of finishing
    the problem, it will suffice.
    """

    component = 0

    for i in range(7):
        for j in range(7):
            if (
                board[i][j].any()
                and (
                    board[i][j + 1].any()
                    or board[i + 1][j].any()
                )
            ) and not (
                board[i][j + 1].any()
                and board[i + 1][j].any()
            ):
                component += 1

    return component


i = 0
while i < 10000:
    game = checkerboard()
    comp_results = test(game)
    comp.append(comp_results)
    i += 1

total = sum(comp)
board_neighbors = 112 * 10000

# while this solution is not exact, it falls with the
# mathematical range mn/8 < E(components) < (m+2)(n+2)/6

print((total / board_spaces * 100), " expected components")
10.93178 expected components

PART 2: Proof of E(c)-range given MxN board:

The total number of neighboring edges given a MxN board is determined by

(m-1)n + (n-1)m

Each square has a probability of 1/2 of being black or white. A random number of shared edges, S, approaches MxN as M and N become sufficiently large.

The expected minimum number of components, E(c) must be

E(c) >= m + n (m-1)(n-1)/8 > mn / 8

This is out lower bound.

With our upper bound twice as large, we simplify:

1/6(m-1)(n-1) + 1/2(n-1) + 4/9(m-1)+1+1/18(m-1)

1/2^(2(n-1)) --> 6E(c) <= (m+2)(n+2)

While not an equation for any given MxN, this is the expected range for any given MxN.