How to make it efficient for group of Gaussian mixture of different size?

The model is

f(\theta_g,X)\sim\sum\limits_{k=1}^{C_g} w_{g,k} N(\mu_{g,k}, \sigma_{g,k}^2);~ g=1:G

where X is data matrix. \theta=[\theta_1,...,\theta_G] is parameter where g denotes group index. \{w_{g,k}, \mu_{g,k}, \sigma_{g,k}^2\}_{k=1}^{C_g} are Gaussian mixture parameters for the g-th group, with C_g Gaussian components.

I tried various ways to code it but it takes 5-15 hours for real data (though only 15 dimensions). Fortunately, it converges well.

How can I make it faster?

Thank you.