ENH - Mixture density class

[AustinRochford]

I often find myself marginalizing over discrete variables and using hacks like

def normal_mixture_rvs(*args, **kwargs):
    w = kwargs['w']
    mu = kwargs['mu']
    tau = kwargs['tau']

    size = kwargs['size']

    component = np.array([np.random.choice(w.shape[1], size=size, p=w_ / w_.sum())
                          for w_ in w])

    return sp.stats.norm.rvs(mu[np.arange(w.shape[0]), component], tau[component]**-0.5)

class NormalMixture(pm.distributions.Continuous):
    def __init__(self, w, mu, tau, *args, **kwargs):
        super(NormalMixture, self).__init__(*args, **kwargs)

        self.w = w
        self.mu = mu
        self.tau = tau

        self.mean = (w * mu).sum()

    def random(self, point=None, size=None, repeat=None):
        w, mu, tau = draw_values([self.w, self.mu, self.tau], point=point)

        return normal_mixture_rvs(w=w, mu=mu, tau=tau, size=size)

    def logp(self, value):
        w = self.w
        mu = self.mu
        tau = self.tau

        return bound(logsumexp(tt.log(w) + (-tau * (value - mu)**2 + tt.log(tau / np.pi / 2.)) / 2.,
                               axis=1).sum(),
                     tau >=0, w >= 0, w <= 1)

to speed up convergence. It seems to me like we could probably add a nice built-in MixtureDensity class to facilitate this flexibly.

Looking for feedback before starting to write this.

[kyleabeauchamp]

Do you recall a good citation that discusses this? Is it in the stan paper?

On Sep 27, 2016 8:05 AM, "Austin Rochford" [email protected] wrote:

I often find myself marginalizing over discrete variables and using hacks
like

def normal_mixture_rvs(_args, *_kwargs):
w = kwargs['w']
mu = kwargs['mu']
tau = kwargs['tau']

size = kwargs['size']

component = np.array([np.random.choice(w.shape[1], size=size, p=w_ / w_.sum())
                      for w_ in w])

return sp.stats.norm.rvs(mu[np.arange(w.shape[0]), component], tau[component]**-0.5)

class NormalMixture(pm.distributions.Continuous):
def init(self, w, mu, tau, _args, *_kwargs):
super(NormalMixture, self).init(_args, *_kwargs)

    self.w = w
    self.mu = mu
    self.tau = tau

    self.mean = (w * mu).sum()

def random(self, point=None, size=None, repeat=None):
    w, mu, tau = draw_values([self.w, self.mu, self.tau], point=point)

    return normal_mixture_rvs(w=w, mu=mu, tau=tau, size=size)

def logp(self, value):
    w = self.w
    mu = self.mu
    tau = self.tau

    return bound(logsumexp(tt.log(w) + (-tau * (value - mu)**2 + tt.log(tau / np.pi / 2.)) / 2.,
                           axis=1).sum(),
                 tau >=0, w >= 0, w <= 1)

to speed up convergence. It seems to me like we could probably add a nice
built-in MixtureDensity class to facilitate this flexibly.

Looking for feedback before starting to write this.

—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
#1401, or mute the thread
https://github.com/notifications/unsubscribe-auth/ABz_q3rU0r6BnBuNFCkinhTEIrbu857hks5quS0sgaJpZM4KHwPR
.

[twiecki]

This would be great, I knew the Stan guys use this to pretty good effect, just didn't know how to do it properly.

[AustinRochford]

@kyleabeauchamp yes, I think it's Chapter 11, Latent Discrete Parameters, of the Stan reference v2.8.0

[AustinRochford]

@twiecki it definitely speeds up mixing drastically for some models. I will start working on this when I get the chance!

[fonnesbeck]

This is essentially what we do in our ZeroInflated* models to eliminate the use of latent indicators. Would be handy to have mixtures of often-used distributions like normals.

[AustinRochford]

@fonnesbeck I was envisioning a base class that could take arbitrary distributions for the mixture components, with subclasses/encapsulating classes for specific common use cases (zero inflated, normal mixtures, etc.).