# Multinomial sampling init val errors out

**URL:** <https://discourse.pymc.io/t/multinomial-sampling-init-val-errors-out/11044>\
**Category:** v5\
**Created:** [December 20, 2022, 6:10am UTC](https://discourse.pymc.io/t/multinomial-sampling-init-val-errors-out/11044 "2022-12-20T06:10:58Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gireesh\_Ramji](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/gireesh_ramji/32/4924_2.png) [@Gireesh\_Ramji](https://discourse.pymc.io/u/Gireesh_Ramji)\
**Post date:** [December 20, 2022, 6:10am UTC](https://discourse.pymc.io/t/multinomial-sampling-init-val-errors-out/11044/1 "2022-12-20T06:10:58Z")

</div>

Hi,

I am trying to use PyMC to estimate HDIs for a couple of variables of interest in the below problem setup.

Let us assume that we have a pool of ids coming from 5 different sources, for which we have access to some aggregated stats. We want to understand which are the sources which are more likely to provide “exclusive” ids, i.e IDs not provided by other sources. More details as below:

1. We get to observe the total number of unique ids in the pool (`n_unique_id_observed`). This includes ids that were contributed exclusively by a source and ids that were contributed by one or more sources.
2. We assume that for each id in the pool, there is a constant probability that this id came exclusively from a given source (one of 5 sources, in this example, `p_exclusive_id_assumption`)
3. We model the number of ids that are from a single source (aka exclusive ids) as a binomial RV (`n_exclusive_id_assumption`)
4. We then assume that for each “exclusive” id, there is a constant probability that this id came from one of the 5 sources (`p_exclusivity_by_source_assumption`)
5. We also get to observe at the source level, the number of ids that came exclusively from each source. We model this as a multinomial RV (`n_exclusive_ids_by_source_assumption`)

The code snippet below generates data as described above:

```auto
n_unique_id_observed = 2500
p_exclusive_id_assumption = 0.10
n_exclusive_id_assumption = np.random.binomial(n=n_unique_id_observed,
                                               p=p_exclusive_id_assumption,
                                               size = 1)
p_exclusivity_by_source_assumption = np.random.dirichlet(np.ones(shape=5))
n_exclusive_ids_by_source_assumption = np.random.multinomial(n_exclusive_id_assumption,
                                                               p_exclusivity_by_source_assumption)

```

Now coming to our real inference challenges:

1. we want to estimate HDIs for `p_exclusivity_by_source`, i.e we are interested in know the range of possibilities which are feasible to answer the question “When we know that an id came from a single source, which source is it most likely to have come from?”

2. we want to estimate HDIs for `p_exclusive_id`, i.e what are the range of possibilities for the fraction of the population which is “exclusive”

To do so, we try to estimate the parameters of a pymc model using the generated data above, as observed. We want to assert that `p_exclusive_id_assumption` and `p_exclusivity_by_source_assumption` lie within the HDI of `p_exclusive_id` and `p_exclusivity_by_source`:

```auto

with pm.Model() as model:
    p_exclusive_id = pm.Beta('p_exclusive_id',
                             alpha=1.0,
                             beta=1.0)
    n_exclusive_id = pm.Binomial('n_exclusive_id',
                             n=n_unique_id_observed,
                             p=p_exclusive_id)
    p_exclusivity_by_source = pm.Dirichlet('p_exclusivity_by_supplier',
                                             a = np.ones(5))
    n_exclusive_ids_by_source = pm.Multinomial('n_exclusive_ids_by_source',
                                             n=n_exclusive_id,
                                             p=p_exclusivity_by_source,
                                             observed = n_exclusive_ids_by_source_assumption)

```

Unfortunately, pm.sample with the above model returns an initval exception:

```auto
SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'p_exclusive_id_logodds__': array(0.), 'n_exclusive_id': array(1158), 'p_exclusivity_by_supplier_simplex__': array([1.60498992, -1.41315054, -0.13382439, 0.18576826])}

Initial evaluation results:
{'p_exclusive_id': -1.39, 'n_exclusive_id': -4.1, 'p_exclusivity_by_supplier': -5.66, 'n_exclusive_ids_by_source': -inf}

```

using the `initval='prior'` keyword arg as suggested on [pm.Multinomial gives Sampling Error - #3 by ricardoV94](https://discourse.pymc.io/t/pm-multinomial-gives-sampling-error/10639/3) does not seem to help either

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [December 20, 2022, 7:02am UTC](https://discourse.pymc.io/t/multinomial-sampling-init-val-errors-out/11044/2 "2022-12-20T07:02:58Z")

</div>

Can you share a reproducible example with data? Glancing at the code I don’t see anything wrong with it, so I suspect a data problem. Either some 0/1 probabilities that don’t match the observations or wrong `n`

---

<div class="post-metadata">

**Author:** ![Gireesh\_Ramji](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/gireesh_ramji/32/4924_2.png) [@Gireesh\_Ramji](https://discourse.pymc.io/u/Gireesh_Ramji)\
**Post date:** [December 20, 2022, 7:47am UTC](https://discourse.pymc.io/t/multinomial-sampling-init-val-errors-out/11044/3 "2022-12-20T07:47:02Z")

</div>

Uploading code to reproduce (same as what I reproduce). Please note that at the moment I am using generated data as `observed` so I doubt it’s a data issue, but maybe I am doing something silly in the way I generate data?

[example.py](https://discourse.pymc.io/uploads/short-url/cltsfUZf7zO7cobZwYJ63o4lQ4r.py) (3.1 KB)

---

<div class="post-metadata">

**Author:** ![Gireesh\_Ramji](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/gireesh_ramji/32/4924_2.png) [@Gireesh\_Ramji](https://discourse.pymc.io/u/Gireesh_Ramji)\
**Post date:** [December 21, 2022, 12:22am UTC](https://discourse.pymc.io/t/multinomial-sampling-init-val-errors-out/11044/4 "2022-12-21T00:22:41Z")

</div>

hi @ricardoV94 does the attached python script help reproduce the error? Thanks so much for looking into this.

---

<div class="post-metadata">

**Author:** ![RavinKumar](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ravinkumar/32/958_2.png) [@RavinKumar](https://discourse.pymc.io/u/RavinKumar)\
**Post date:** [January 18, 2023, 2:42pm UTC](https://discourse.pymc.io/t/multinomial-sampling-init-val-errors-out/11044/5 "2023-01-18T14:42:45Z")

</div>

Hey @Gireesh_Ramji  
To better answer this question here’s some things things that will help. One is more cersion numbers for the major packages, at least PyMC and Pytensor.

For debugging here’s two strategies. Have you tried fixing certain distributions to specific values and seeing if sampling workings then? Perhaps `n_exclusive_id_assumption`

Another is passing in specific test point values and seeing if you get a logp and gradient back in regions you expect. This is a screenshot from my book. The API is PyMC3, so I don’t think the code 1:1 will work, but the idea is the same. Check if the model is returning reasonable values in places you think it should.

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/9/975b36172bd0e2f3bd25cc7c151fa098c00e50d2.png)

---

<div class="post-metadata">

**Author:** ![Gireesh\_Ramji](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/gireesh_ramji/32/4924_2.png) [@Gireesh\_Ramji](https://discourse.pymc.io/u/Gireesh_Ramji)\
**Post date:** [January 18, 2023, 10:56pm UTC](https://discourse.pymc.io/t/multinomial-sampling-init-val-errors-out/11044/6 "2023-01-18T22:56:47Z")

</div>

Thank you @RavinKumar - the versions are:  
PyMC: 4.2.2  
Aesara: 2.8.7

So I did find a fix and that is the observation that `n_exclusive_id` may also be constrained by providing its observed value as `n_exclusive_id_by_source_assumption.sum()`

So essential replacing

```auto
    n_exclusive_id = pm.Binomial('n_exclusive_id',
                             n=n_unique_id_observed,
                             p=p_exclusive_id)

```

with

```auto
    n_exclusive_id = pm.Binomial('n_exclusive_id',
                             n=n_unique_id_observed,
                             p=p_exclusive_id,
                             observed = n_exclusive_id_by_source_assumption.sum())

```

I guess I was expecting that providing this additional constraint to the model was not necessary as I thought the model would have converged to this value of n\_exclusive\_id\_by\_source\_assumption.sum() for n\_exclusive\_id, on account of the specification that n\_exclusive\_id is the number of trials for n\_exclusive\_id\_by\_source.

However it appears that the initial value checks do not have any way to infer this implicit constraint and hence is unable to reconcile the initial (internally proposed) value for `n_exclusive_id` with the observed value provided for `n_exclusive_id_by_source`.

Been getting tripped up quite frequently with initial value checks - would really appreciate being pointed to a resource that will help me understand the overall motivation and desired impact of running the initial value checks.