Extended Monty Hall problem

Questions
1

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?

thanks.

2

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

2 Likes
3

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