1
$\begingroup$

I'm sure this is a simple problem, but I'm not quite sure if I'm doing this right.

Let me preface this with I am a computer-science person, and will need a tractable solution that I can program.

In most programming languages, you have a function rand() that will return a real number from 0 to 1. This is a uniform distribution.

However, I am trying to create a galaxy of stars. Thus they will be more dense toward the center, and less dense in the outer reaches of the galaxy.

I also want to set the maximum densities toward the center and minimum densities toward the outside edge. So for example, 5% density at the outer edge, and 60% density at the inner edge. And in between would be some scale of the cdf of the normal. (From which I will get and use the equation from wiki here.)

So I would do this = f(rand(), edgeDensity, coreDensity)

I was thinking of taking $coreDensity - edgeDensity$ and adjust the scale of the normal distribution to be

$f(x, e, c) = [\phi(x) \cdot (c-e)] + e$,

where we have the following:

  • $0 \leq e < c \leq 1$

  • $0 \leq x \leq 1$

  • $\displaystyle \phi(x) = \frac{1}{2}\left[1+\mathrm{erf}\left(\frac{x-\mu}{(2\sigma^2)^{1/2}}\right)\right]$

But this is where I get stuck. I'm not sure which values I should use for $\sigma$ and $\mu$ as well as I know $\mathrm{erf}$ is some error function, but is there a simple non integral solution to it?

Or is there a better way to do this?

  • 1
    Nope, there is no simple closed form. Liouville proved this in the early-to-mid 1800's. There are very good numerical approximations, though. For starting points, see [here](http://stats.stackexchange.com/questions/7200/evaluate-definite-interval-of-normal-distribution) and [here](http://stats.stackexchange.com/questions/9501/is-it-possible-to-analytically-integrate-x-multiplied-by-the-lognormal-probabil), and, **especially**, W. J. Cody's paper **[here](http://www.ams.org/journals/mcom/1969-23-107/S0025-5718-1969-0247736-4/S0025-5718-1969-0247736-4.pdf)**.2011-05-03
  • 0
    what is the $R_{lm}()$ function in Cody's paper.2011-05-03
  • 0
    A rational function of degree $\ell$ and $m$ in the numerator and denominator, respectively, as stated at the bottom of the first page. Coefficients for the rational functions are given on the subsequent pages.2011-05-04
  • 0
    So it looks like the $$R_{lm}() = -100 \log_{10} max\Bigl(\dfrac{f(x) - f_{lm}(x)}{f(x)}\Bigr)$$ but it says that $f(x) = erf(x)$ which would imply a recursive element? or am I reading this wrong?2011-05-04
  • 0
    The error function $\mathrm{erf}(x)$ certainly cannot be expressed in terms of "elementary" functions (but there are [ways](http://math.stackexchange.com/questions/97) to compute it on a computer); on the other hand, you seem to want to have a method for generating normally-distributed (pseudo)random numbers, in which case, there is a way to do it without the $\mathrm{erf}$ rigamarole. If this is what you want, I'll write something up.2011-05-04
  • 0
    That would be fantastic. I was going down this path, because it seemed most reasonable to me, but if there is a better one, I'd take that for sure.2011-05-04
  • 1
    For generating normal random variates: **[Box-Muller transform](http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform)**. It's not the fastest algorithm out there (cf. [ziggurat](http://en.wikipedia.org/wiki/Ziggurat_algorithm)), but it's still quite fast and very easy to implement.2011-05-04

0 Answers 0