Extended Monty Hall problem

As I was going through the pymc 4.0 tutorial, I thought I could extend the Monty Hall problem to having multiple doors and really observe the addition of the probabilities.
I started by defineing this helper function that report array of probabilities of left over doors at all scenario:

def left_over_door(num_of_doors: int = doors, chosen: int = 0) -> np.ndarray:
    res = np.zeros(shape=(num_of_doors, num_of_doors))
    np.fill_diagonal(res, val=1)
    res[chosen] = np.full(shape=(1, num_of_doors), fill_value = 1./(num_of_doors-1))
    res[chosen, chosen] = 0
    return res

then I proceed to the model inference, by using at.take instead of at.switch in the example.

doors = 3
with pm.Model() as monty_multiple_doors_model:
    # prior setup
    prize = pm.DiscreteUniform("prize", 0, doors-1, initval = 0)
    p_not_open = at.take(left_over_door(doors, chosen=0), indices=prize, axis=0)
    # likelihood
    not_opened = pm.Categorical('not_opened', p_not_open, observed=2)

however, the model failed at sampling step. Here’s the error:

IndexError                                Traceback (most recent call last)
IndexError: index out of bounds
Apply node that caused the error: Subtensor{int64}(TensorConstant{[[0.  0.5 .. 0.  1. ]]}, ScalarFromTensor.0)
Toposort index: 8
Inputs types: [TensorType(float64, matrix), Scalar(int64)]
Inputs shapes: [(3, 3), ()]
Inputs strides: [(24, 8), ()]
Inputs values: ['not shown', 4]
Outputs clients: [[Elemwise{ge,no_inplace}(Subtensor{int64}.0, TensorConstant{(1,) of 0}), Sum{acc_dtype=float64}(Subtensor{int64}.0), Elemwise{true_div,no_inplace}(Subtensor{int64}.0, InplaceDimShuffle{x}.0)]]
RuntimeError: Chain 0 failed.

I find it stange, 'cause the way I specify my model shouldn’t lead to any index error. Could it be the way I call the aesara API which lead to this? If it is, how should I fix/war my code to make the simulation work?


The reason it failed during sample is that, PyMC assign Metropolis as sampler to the random variable prize, which can generate proposal out of the range [0, 2], which cause the index error.
Changing prize to

    prize = pm.Categorical("prize", np.ones(3)/3)

should work


after implementing the solution, the sampler changed to CategoricalGibbsMetropolis. I guess, from what the name suggests, is a sampler that out-of-range proposals are never proposed.
and from the API, I think that at.switch (used in the example) can accept out-of-range proposals while the at.take I’ve used does not. so I think this concludes my problem.

1 Like