Prevent infinite loop when a prior is discrete AND using Multivariate Normal Transition

[Edderic]

Hi!

First of all, I'd like to thank everyone who's contributed to this project. I'm excited to use it for research!

I'd like to get your thoughts about the issue I discovered. Over the last few days, I had spent a lot of time trying to figure out why my code would never progress past the second population even though the population size was set to 10 (also making use of default settings). Good news is that I think I finally figured it out. The scenario I'm running into is that I'm getting stuck in an infinite loop because there's a mismatch between using discrete priors (e.g. RV('randint', 0, 10)) and using the MultivariateNormalTransition, which assigns continuous values instead of discrete, leading to an infinite loop (of bad proposals).

In the _generate_valid_proposal method, we get samples from MultivariateNormalTransition by default:

pyABC/pyabc/smc.py

Line 600 in f6941cf

theta_ss = transitions[m_ss].rvs()

We then figure out if the model and parameter priors are plausible. We only get out of the while loop if the model and parameter priors have a non-zero probability:

if (model_prior.pmf(m_ss)
      * parameter_priors[m_ss].pdf(theta_ss) > 0):
       return m_ss, theta_ss

Let's look at the parameter_priors[m_ss].pdf(theta_ss) part. See this link to that part of the codebase:

            for key, val in x.items():
                try:
                    # works for continuous variables
                    res.append(self[key].pdf(val))
                except AttributeError:
                    # discrete variables do not have a pdf but a pmf
                    res.append(self[key].pmf(val))
            return reduce(lambda s, t: s * t, res)

When self[key] is discrete, we get to the except section. So we'll evaluate self[key].pmf(val). The problem, though, is that using MultivariateNormalTransition gives a continuous value, so doing self[key].pmf(val) gives us 0. return reduce(lambda s, t: s * t, res) would then be guaranteed to be 0, so we'll be stuck in an infinite-loop.

I recommend changing res.append(self[key].pmf(val)) to res.append(self[key].pmf(int(round(val)))) so that we get useful values. Another option is to instead raise an error to give the user an idea that using MultivariateNormalTransition doesn't make sense when discrete priors are used. Do you think there's a better way that I haven't discussed? Would love to hear your thoughts on this. Please let me know if you need any more details. Thank you!

Regards,
Edderic

[neuralyzer]

Hi. I think the rounding would not work for all discrete variables. Disctete variables are not necessarily integer. We might really need another transition here? For example model selection uses another trandition already. Also, integers and categorical variables might both be encoded as integers. The kind of transition they require is however likely different. ⁣Sent from BlueMail ​

…

On Mar 24, 2019, 04:51, at 04:51, Edderic Ugaddan ***@***.***> wrote: Hi! First of all, I'd like to thank everyone who's contributed to this project. I'm excited to use it for research! I'd like to get your thoughts about the issue I discovered. Over the last few days, I had spent a lot of time trying to figure out why my code would never progress past the second population even though the population size was set to 10 and assuming all other default settings. Good news is that I think I finally figured it out. The scenario I'm running into is that I'm getting stuck in an infinite loop because there's a mismatch between using discrete priors (e.g. `RV('randint', 0, 10)`) and using the `MultivariateNormalTransition`, which assigns continuous values instead of discrete, leading to an infinite loop (of bad proposals). In the `_generate_valid_proposal` method, we get samples from `MultivariateNormalTransition` by default: https://github.com/ICB-DCM/pyABC/blob/f6941cfe347d761e2ab3f1463678e28da4a173f4/pyabc/smc.py#L600 We then figure out if the model and parameter priors are plausible. We only get out of the `while` loop if the model and parameter priors have a non-zero probability: ```python if (model_prior.pmf(m_ss) * parameter_priors[m_ss].pdf(theta_ss) > 0): return m_ss, theta_ss ``` Let's look at the `parameter_priors[m_ss].pdf(theta_ss)` part. [See this link to that part of the codebase](https://github.com/ICB-DCM/pyABC/blob/9917eed639c3f60d64bae9dafa4e04e5bb922d83/pyabc/random_variables.py#L449): ```python for key, val in x.items(): try: # works for continuous variables res.append(self[key].pdf(val)) except AttributeError: # discrete variables do not have a pdf but a pmf res.append(self[key].pmf(val)) return reduce(lambda s, t: s * t, res) ``` When `self[key]` is discrete, we get to the `except` section. So we'll evaluate `self[key].pmf(val)`. The problem, though, is that using `MultivariateNormalTransition` gives a continuous value, so doing `self[key].pmf(val)` gives us 0. `return reduce(lambda s, t: s * t, res)` would then be guaranteed to be 0, so we'll be stuck in an infinite-loop. I recommend changing `res.append(self[key].pmf(val))` to `res.append(self[key].pmf(int(round(val))))` so that we get useful values. Another option is to instead raise an error to give the user an idea that using `MultivariateNormalTransition` doesn't make sense when discrete priors are used. Do you think there's a better way that I haven't discussed? Would love to hear your thoughts on this. Please let me know if you need any more details. Thank you! Regards, Edderic -- You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub: #75

[Edderic]

Hi @neuralyzer,

Thanks for responding so quickly. I think your point definitely makes sense. For categorical ordinal variables, the multivariate normal transition with my proposed rounding might still work. However, in cases where it's categorical but not ordinal, it doesn't make sense at all.

We might really need another transition here? For example model selection uses another transition already.

You're talking about this method right?

def _get_discrete_rv(self, m):
        p_stay = self.probability_to_stay
        p_move = (1 - p_stay) / (self.nr_of_models - 1)
        probabilities = [p_stay if n == m else p_move
                         for n in range(self.nr_of_models)]
        return RV('rv_discrete', values=(range(len(probabilities)), probabilities))

So when specifying a prior for discrete variables, does it make sense to require an explicit argument for the ordinality? (e.g. RV('randint', 0, 10, ordinal=True). Then, when generating a proposal, we make use of the ordinal field to decide which transition to use when the transition isn't explicitly specified? Do you think a transition function similar to above would work?

[yannikschaelte]

Hi, since you use a discrete prior, I assume your parameter space is indeed discrete? In that case, the MultivariateNormalTransition might actually not be the best choice. If the rounding works for you (you can simply subclass the prior Distribution to apply the rounding), then that should be fine though.

I had at some point implemented a discrete transition based on random walks, assuming integer ordinal parameters, but so far we have not found use for it yet. I uploaded it here now, in the discrete_transition branch. It is certainly not as tuned as the Gaussian transitions, but it might be worth a try.

Agreed that there should be some sort of warning when using a discrete distribution and continuous transitions. We will think about how to best to detect and do that, thanks!

[Edderic]

Hi, since you use a discrete prior, I assume your parameter space is indeed discrete? In that case, the MultivariateNormalTransition might actually not be the best choice. If the rounding works for you (you can simply subclass the prior Distribution to apply the rounding), then that should be fine though.

Yes that makes sense! For the project I'm working on, the variables I've thought about so far are ordinal, so I could technically just simplify the problem and still have some good-enough approximate answer by switching to continuous variables. However, I think subclassing the prior Distribution might be useful, even in other cases. For example, one might have a constant prior. In that case, we don't really want to transition anywhere else, so for the .rvs method, we should just return the constant.

Agreed that there should be some sort of warning when using a discrete distribution and continuous transitions. We will think about how to best to detect and do that, thanks!

Yeah that would be awesome! It certainly would have saved me weeks of wondering if the infinite-loop issue was in my model and not in PyABC 😄

[Edderic]

I had at some point implemented a discrete transition based on random walks, assuming integer ordinal parameters, but so far we have not found use for it yet. I uploaded it here now

@yannikschaelte I tried clicking on the link you gave but it shows me "Page not found"?

[yannikschaelte]

Oh yes, sorry. That is because that branch got deleted. It is in master branch now: https://github.com/ICB-DCM/pyABC/blob/master/pyabc/transition/randomwalk.py

(To clarify: Since we so far unfortunately never used a discrete distribution, the implementation is rather naive, and most certainly can be improved.)

[Edderic]

Thanks, @yannikschaelte!

In the case where one has a mix of continuous and discrete priors, we'll have to have more than one transition function, such as the MultivariateNormalTransition and the RandomWalk, right? So the ABCSMC class will have to deal with that change? (e.g. expect maybe a ModelTransition object that's a wrapper for a chosen continuous and discrete transition?)

[yannikschaelte]

Hi! Yes, when there is the necessity to distinguish in the treatment between dimensions (e.g. some discrete, some continuous), one would, just as for the prior Distribution, need to write an own Transition which is then passed to ABCSMC. This new "AggregatedTransition" could hold e.g. a MultivariateNormalTransition and a DiscreteTransition as attributes, and pass the reduced input on to each of them for the respective coordinates, and then aggregate the result (e.g. to get a new proposal step, or multiply the probabilities, assuming independence, to get the proposal probability).
There might be a theoretical problem though with how to define the probability for a mixed discrete/continuous proposal probability, not sure about that.

I think writing such an AggregatedTransition should be pretty straightforward for you. I will make an issue on having one inside pyabc too to offer a generic implementation if possible at some point.

[yannikschaelte]

Closed now due to inactivity. Separate issues were created for the discussed points.

Please feel free to re-open if there are questions left :) .