Dear all,

I am currently working with Finite Mixture Models, which, as you may know, involve a high number of parameters. An MCMC estimation of these parameters may take a long time and I have been recommended to use a Hamiltonian MC sampler.

I have just started reading about RStan and how to define my model. My Mixture Model does not involve conventional distributions so I have to declare it in the functions block. My probability distribution function is the Einasto profile (this function is used in astrophysics for modeling galaxies and halos):

exp(-dn*((dist/rs)**(1/n)-1))

where,

dist = distances of the observations y to the center of the halo (mu)

mu = center of the halo (parameter of the distribution)

rs = radius of the halo (parameter of the distribution)

n = shape of the halo (parameter of the distribution)

So far, so good, but the problem comes now. The quantity dn is defined as the quantity satisfying the equation

gamma(3*n) = 2*gammainc(3*n,dn)

where n is the mentioned shape of the distribution (parameter). gamma and gammainc are the Gamma and the Incomplete Gamma functions respectively.

When I wrote the code in R for the first time I had to solve the equation every time I evaluated the Einasto profile using an optimization routine. How can I do this with RStan? Can I call an optimization routine inside the functions block? Is there any?

Many thanks in advance for your input.