I'm have to implement Gaussian Normal distribution in a programming language and now I'm trying to understand what is a Gaussian Normal distribution.

I have read its Wikipedia entry, but sorry, there are a lot of formula and a lot of things that I don't understand.

If I want to calculate a Gaussian distribution with zero mean and standard deviation $\sigma$, or $N(0,\sigma)$, do I need to implement the probability density function? $$p(\Delta x_i) = \frac{1}{\sqrt{2\pi\sigma}}e^{-(\Delta x_i)^2/2\sigma^2}$$

I have copied the above formula from the book Introduction to Evolutionary Computing which is the one I'm studying.

In other words: I need to calculate $x_i^{\prime} = x_i + N(0, \sigma)$ but I don't know how to calculate $N(0, \sigma)$.

I'm going to use C# to implement it and I don't have any other resources to do it, just Microsoft's .Net Framework.

  • 1
    $\begingroup$ You won't really need the density, just the CDF. Try reading an introductory statistics text. $\endgroup$ Jan 20, 2018 at 16:07
  • $\begingroup$ Well, that is the formula that appears in the book, this is why I'm asking. If they wanted me to use the CDF formula, I think they would have put it. $\endgroup$
    – VansFannel
    Jan 20, 2018 at 16:11
  • 1
    $\begingroup$ A random variable $X$ having a Gaussian distribution with mean zero and sd $\sigma$, usually denoted by $X\sim\mathcal N(0,\sigma^2)$, has the density function $f(x)=\frac{1}{\sqrt{2\pi}\sigma}e^{-(x^2/2\sigma^2)}$ where $x\in\mathbb{R}$. $\endgroup$ Jan 20, 2018 at 16:15
  • $\begingroup$ @StubbornAtom thanks but I need to calculate $N(0,\sigma)$ and I don't know how to do it. $\endgroup$
    – VansFannel
    Jan 20, 2018 at 16:23
  • 1
    $\begingroup$ It's not quite as simple as implementing the PDF. I can share my code in C for doing this later. For now the following link should be of use stackoverflow.com/questions/218060/random-gaussian-variables . $\endgroup$
    – Spencer
    Jan 20, 2018 at 16:38

1 Answer 1


First I will sketch out a quick "Physicists Proof" of the method. This approach was taught to me by the Astrophysicist James Fry.

To generate samples of a random variable $X$ which has a cumulative distribution function $F(x)$ you apply the inverse of $F$ to samples from a uniform random variable on the interval from $[0,1]$.

This is simple enough, but the cumulative distribution function for a normal random variable doesn't have a simple closed form. If we want to stick with elementary functions, which we are likely to have well implemented in our libraries, we should avoid directly computing this inverse. What we do instead is re-express the problem in terms of a pair of random variables whose distribution functions are more easily inverted in terms of elementary functions.

Let $X$ and $Y$ be independent normal random variables which obey $N(0,\sigma)$. These have a cumulative distribution function,

$$ F(x,y) = \int_{-\infty}^{x} \int_{-\infty}^{y} dx' dy' \ \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{x'^2}{2\sigma^2} } \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{y'^2}{2\sigma^2} } $$

now we are going to change variables to polar coordinates, but I am not going to pay too much attention to the specific bounds. Recall that $x=r\cos(\theta)$ and $y=r\sin(\theta)$ and $dx\ dy = rdr\ d\theta$

$$ F(r,\theta) = \int_{0}^{r} \int_{0}^{\theta} \frac{r'dr' d\theta'}{2\pi\sigma^2} \ e^{-\frac{r'^2}{2\sigma^2} }$$

this integral can be factored,

$$ F(r,\theta) = \int_{0}^{r} \frac{r'dr' d\theta'}{\sigma^2} \ e^{-\frac{r'^2}{2\sigma^2} }\int_{0}^{\theta} \frac{d\theta'}{2\pi}$$

the $r$ integral can be modified by a change of variables $t=r^2/2\sigma^2$.

$$ F(t,\theta) = \int_{0}^{t} dt' \ e^{-t'} \int_{0}^{\theta} \frac{d\theta'}{2\pi}$$

Each integral represents the cumulative distribution function of a random variable ($T$ and $\Theta$). These distribution functions are easy to compute.

$$F_T(t) = 1-e^{-t} \qquad F_\Theta(\theta) = \frac{\theta}{2\pi}$$

Each of these functions can easily be inverted,

$$F_T^{-1}(t) = -\ln(1-t) \qquad F_\Theta^{-1}(\theta) = 2\pi\theta,$$

so we can easily generate samples for the random variables $T$ and $\Theta$. These samples are then used to compute $r$ and subsequently $x$ and $y$.

The procedure is as follows.

1) By some means, generate a pair of random numbers which obey an uniform distribution on the interval [0,1]. Call these numbers $u_1$ and $u_2$.

2) Compute numbers $t$ and $\theta$ by, $t=-\ln(1-u_1)$ and $\theta=2\pi u_2$.

3) Compute a number $r$ by, $r=\sqrt{2\pi\sigma^2 t} $.

4) Compute a pair of numbers $x$ and $y$ by, $x=r\cos(\theta)$ and $y=r\sin(\theta)$.

5) $x$ and $y$ are samples from a $N(0,\sigma)$ distribution. If you want samples from a $N(\mu,\sigma)$ distribution then add $\mu$ to each of $x$ and $y$.

Below is code I have in C for generating these numbers. To compile this code with a GNU C compile simply type "gcc name_of_file_with_code.c -lm".

#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <complex.h>

double rUniform(double, double);
double complex rNormal(double, double);

const int ITMAX=1000;

int main(){

  for( int i=0; i < ITMAX; i++){
    complex double z = rNormal(5,1); 
    printf("%f\t%f\n", creal(z),cimag(z));

double rUniform(double a, double b){
  return rand()*(b-a)/RAND_MAX+a;

double complex rNormal(double mu, double sigma){
  // Generate a pair of random numbers which are
  // sampled from a uniform distribution on the
  // interval [0,1]. 
  double u1 = rUniform(0,1);
  double u2 = rUniform(0,1); 

  double pi = acos(-1); // Lazy way to get pi=3.14... from arccos(-1)

  // Generte samples for the variables Theta and T by
  // applying the inverses of the cumulative distribution
  // functions to uniform random variables. 
  double theta = 2*pi*u1;
  double t = -log(1-u2);

  // Generate r from t, recall that t=r^2/(2*sigma^2)
  // so we can get r from sqrt(2*sigma^2*t)
  double r = sqrt(2*sigma*sigma*t);

  // Generate the cartesian coordinates x and y from the polar
  // coordinates r and theta. Recall that these x and y are the
  // variables which actually obey N(mu,sigma).
  double x = r*cos(theta);
  double y = r*sin(theta);

  // We can return the pair in a complex variable z.
  double complex z =  (x+mu)+I*(y+mu);

  return z;  
  • $\begingroup$ Thanks a lot. So, If I have understood what is a normal distribution, given a mean, in my case with value 0, and a standard deviation, $N(0,\sigma)$ will return random numbers from the Gaussian distribution, isn't it? $\endgroup$
    – VansFannel
    Jan 21, 2018 at 6:52
  • $\begingroup$ I'm a little confused by your question. Generally Gaussian and Normal are used interchangably to describe the same distribution. $\endgroup$
    – Spencer
    Jan 21, 2018 at 15:13
  • $\begingroup$ Yes, I'm talking about normal distribution. If I run $N(0,\sigma)$ three times it will return three random numbers, isn't it? $\endgroup$
    – VansFannel
    Jan 21, 2018 at 19:45
  • $\begingroup$ The notation $N(0,\sigma)$ is the name of the distribution. The function in my code rNormal, will return random numbers sampled from this distribution with the appropriate mean and standard deviation. So your understanding is basically correct. $\endgroup$
    – Spencer
    Jan 21, 2018 at 21:12
  • $\begingroup$ I have seen your profile and I would like to study mathematics. What do you recommend me? A book? A MOOC? Thanks. $\endgroup$
    – VansFannel
    Jan 22, 2018 at 16:44

Not the answer you're looking for? Browse other questions tagged .