---
abstract: |
    Algorithm and code to produce sequences whose members obey Gaussian distribution function is reported. Discreet and limited number of groups are defined in the distribution function, where each group is represented only with one value instead of a range of value. The produced sequences are also checked back whether they still fit the discreet distribution function. Increasing of number of particles $N$ increases the value of correlation coefficient $R^2$, but increasing number of groups $M$ reduces it. Value $R^2 = 1$ can be found for $N = 1000000$ at least with $M = 5000$ and for $M = 10$ at least with $N = 1000$.\
    [**Keywords:**]{} gaussian distribution, random sequence, algorithm, code.
author:
- |
    Sparisoma Viridi\*\
    Nuclear Physics and Biophysics Research Division\
    Institut Teknologi Bandung, Bandung 40132, Indonesia\
- |
    Veinardi Suendo\
    Physical Chemistry Research Division\
    Institut Teknologi Bandung, Bandung 40132, Indonesia\
- '\*dudung@fi.itb.ac.id'
bibliography:
- 'manuscript.bib'
title: 'How to produce discreet Gaussian sequences: Algorithm and code'
---

Introduction
============

Gaussian distribution function plays important role in many fields of science, such as in mathematical modeling [@Gershenfeld_2002], in physical sciences [@Boas_1983], in quantum chemistry [@Atkins_1970], with integral in nuclear physics [@Duderstadt_1976], and in semiconductor devices [@Sze_1985]. Then a need comes up how a sequence, that its members obey Gaussion distribution function, could be produced, since it is needed, for example in molecular dynamics simulatons [@Sustini_2011]. A procedure to produce the sequences is presented in algorithm and C++ code.

Gaussian distribution function
==============================

Gaussian or normal distribution function can be represented in the form of

$$\label{eq01}
f(z) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left[-\frac{(z - \mu)^2}{2\sigma^2}\right]$$

where $\mu$ is the average of $z$ and $\sigma$ is the width of normal distribution curve. The factor in front of right side of Equation (\[eq01\]) is due to normalization of $f(z)$ integral

$$\label{eq02}
\int_{-\infty}^{\infty} f(z) dz = 1.$$

Variable $z$ is a certain parameter that obeys Gaussian distribution function, it can be particle velocity, particle diameter, or particle mass.

Proof of normalization
----------------------

Equation (\[eq02\]) can be proved using

$$\label{eq03}
\int_{-\infty}^{\infty} \exp{(-x^2)} dx = \sqrt{\pi},$$

so that

$$\begin{aligned}
\nonumber
\int_{-\infty}^{\infty} f(z) dz = \int_{-\infty}^{\infty} \frac{1}{\sigma\sqrt{2\pi}} \exp\left[-\frac{(z - \mu)^2}{2\sigma^2}\right] dz \\
\nonumber
= \frac{1}{\sigma\sqrt{2\pi}} \int_{-\infty}^{\infty} \sigma\sqrt{2} \exp\left[-\frac{(z - \mu)^2}{2\sigma^2}\right] d\left(\frac{z - \mu}{\sigma\sqrt{2}}\right) \\
\nonumber
= \left(\frac{1}{\sigma\sqrt{2\pi}}\right)\left(\sigma\sqrt{2}\right) \int_{-\infty}^{\infty} \exp\left[-\frac{(z - \mu)^2}{2\sigma^2}\right] d\left(\frac{z - \mu}{\sigma\sqrt{2}}\right) \\
\nonumber
= \left(\frac{1}{\sigma\sqrt{2\pi}}\right)\left(\sigma\sqrt{2}\right) \left(\sqrt{\pi}\right) = 1.\end{aligned}$$

Meaning of $\mu$ and $\sigma$
-----------------------------

Peak of $f(z)$ is located at $z = \mu$ with value

$$\label{eq04}
f_{\max}(z) = f(\mu) = \frac{1}{\sigma\sqrt{2\pi}}$$

and at $z = \mu \pm \frac12 \sigma$ it gives

$$\label{eq05}
f\left(\mu \pm \frac12 \sigma\right) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left[-\frac{1}{8}\right] \approx \frac{0.8825}{\sigma\sqrt{2\pi}}.$$

An example of $f(z)$ is given in Figure \[fg01\].

![\[fg01\] An Gaussian distribution function with $\mu = 0.5$ and $\sigma = 0.25 / \sqrt{2\pi}$.](fig01.eps){width="9cm"}

Number of particles
-------------------

Suppose that there are $N$ particles in a system, that number of particles $N(z)$ who has property of $z$ is defined by

$$\label{eq06}
N(z) = \frac{N}{\sigma\sqrt{2\pi}} \exp\left[-\frac{(z - \mu)^2}{2\sigma^2}\right],$$

where according to Equation (\[eq02\]) it must hold that

$$\label{eq07}
\int_0^{\infty} N(z) dz = N.$$

In this case, it is considered that property $z$ has only positive value.

Discretization of distribution function
=======================================

It is imposible even with nowadays most advanced computer facilities to produce continue number of particles in order of one mole, which equals to about $10^{23}$ particles. In this report only small number of particles is considered. The distribution function is also simplified by dividing it into limited and discreet groups of particles. Within each group there is only one value (certain property of particle) which represents the group instead of a range of value from minimum to maximum value of the group.

Discreet groups
---------------

Suppose that there is $M$ groups of particles with equal width $\Delta z$, which group the total number of particles $N$. First step is how to find $z_{\min}$ and $z_{\max}$ where at these values $N(z)$ can be considered zero. Since we deal with particles than it is more simple to use the [int()]{} function which returns the integer value of $N(z)$. It means that $f(z)$ is considered zero when $N(z) = 1 - \epsilon$, then

$$\begin{aligned}
\label{eq08}
N(z) = 1 - \epsilon, ~z < \mu \Rightarrow z = z_{\rm \min}, \\
\label{eq09}
N(z) = 1 - \epsilon, ~z > \mu \Rightarrow z = z_{\rm \max},\end{aligned}$$

with $\epsilon$ a small defined value. Then width $\Delta z$ can be found through

$$\label{eq10}
\Delta z = \frac{z_{\rm \max} - z_{\rm \min}}{M}.$$

Group $i$ is represented by $z_i$, which is

$$\label{eq11}
z_i = z_{\rm \min} + \left(i - \frac12 \right)\Delta z, ~i = 1, 2, .., M - 1, M.$$

Member of each group
--------------------

As it has been declared previously, in group $i$ there is only one value of $z$ which is $z_i$. It is only for the sake of simplicity. Each group has number of particles that must obey the Gaussian distribution function. Number of particle in each group is

$$\label{eq12}
N_i = \left(\frac{N}{N'}\right) {\rm int}[N(z_i)].$$

Since there is a round down process (throug the [int()]{} function) for each group in order to find $N_i$ from $N(z_i)$ then it can be concluded that

$$\label{eq13}
\sum_{i = 1}^M N_i \le N,$$

a difference that deviates the discreet groups of particles from the Gaussian distribution function. The factor in front of right side of Equation (\[eq12\]) is due to discreet number of particle groups.

Algorithm to group the particles
--------------------------------

An algorithm of implementation of Equation (\[eq08\]) - (\[eq12\]) can be summerized as follow

1.  [\[st11\]]{} start

2.  [\[st12\]]{} determine mu and sigma for distribution function N(z)

3.  [\[st13\]]{} determine epsilon

4.  [\[st14\]]{} set z = mu

5.  [\[st15\]]{} using root finding algoritm find root of N(z) - (1 - epsilon) = 0 in range z &lt; mu, it is named as zmin

6.  [\[st16\]]{} set z = mu

7.  [\[st17\]]{} using root finding algoritm find root of N(z) - (1 - epsilon) = 0 in range z &gt; mu, it is named as zmax

8.  [\[st18\]]{} determine number of group M

9.  [\[st19\]]{} calculate group width dz using Equation (\[eq10\])

10. [\[st1A\]]{} determine zi using Equation (\[eq11\]) for all M groups

11. [\[st1B\]]{} deterimine number of group i using Equation (\[eq12\])

12. [\[st1C\]]{} calculate N’ and normalize Ni with it

13. [\[st1D\]]{} stop

The sequences
=============

In group $i$ there are $N_i$ particles which has a property $z_i$. The property can be velocity, mass, diameter, charge, or other physical properties. And there are $M$ groups of particles. It means, when all the particles are lined in order to make sequences there will be $S$ ways to rearrange the particles order. If the particles are distinguishable

$$\label{eq14}
S_{\rm distinguishable} = N!$$

and when there are indistinguishable

$$\label{eq15}
S_{\rm indistinguishable} = \frac{N!}{\prod_{i = 1}^M N_i!}.$$

The later means that particles at the same group are identical, which means the particles are identify only by their property $z_i$.

The zeroth sequence
-------------------

The easiest way to buid the sequence is by lining the particle from each group in incremental order, such as

$$\label{eq16}
z_1, z_1, z_2, z_2, z_2, z_2, z_3, z_3, .., z_M, z_M.$$

This sequence is named as the zeroth sequence.

Other sequences
---------------

The sequences beside zeroth sequence can be generated by permutating zeroth sequence. Number of sequences can be produced is according to Equation (\[eq14\]) and (\[eq15\]). In this report we propose a mechanism to generate a sequence from zeroth sequence by using [random()]{} and [swap()]{} function which is already built-in in C++. The algorithm is as follow

1.  [\[st21\]]{} start

2.  [\[st22\]]{} determine seed for random generator

3.  [\[st23\]]{} set the generator with the seed

4.  [\[st24\]]{} get the zeroth sequence that contains N particles

5.  [\[st25\]]{} particle number i = 1

6.  [\[st26\]]{} generate an integer number between 1 and N, say j

7.  [\[st27\]]{} swap value of particle i and j

8.  [\[st28\]]{} increase value of i by 1

9.  [\[st29\]]{} if i still less than or equal to M go to Step \[st26\]

10. [\[st2A\]]{} stop

Since random number generated by C++ random generator depends on the seed, than the sequence is reproducible. It means that the seed is as an idenfier to the sequence.

Error
=====

The sum of generated value of $N_i$ for each group $i$ in a sequence will be less than total number of particle as given by Equation (\[eq13\]), which means an error. This error can be calculated using a common correlation coefficient $R^2$ formulation

$$\label{eq17}
R^2 = 1 - \frac{SS_{\rm err}}{SS_{\rm tot}},$$

where

$$\begin{aligned}
\label{eq18}
SS_{\rm err} = \sum_i [N_i - N(z_i)]^2, \\
\label{eq19}
SS_{\rm tot} = \sum_i (N_i - \overline{N}_i)^2, \\
\label{eq20}
\overline{N}_i = \frac{1}{N'} \sum_i N_i,\end{aligned}$$

with $N'$ is total number of generated particles

$$\label{eq21}
N' = \sum_{i = 1}^M N_i.$$

Equation (\[eq17\]) - (\[eq21\]) will be used the next section to calculate the error in produced sequences.

Results and discusion
=====================

An illustration for two discreet Gaussian distribution function is given in Figure \[fg02\], which is produced by our program [gaussg]{}. It has been found that the value $N'$ shown in Equation (\[eq12\]) can not be used in the continue function to fit the discreet values. Then, the new fitting function will be

$$\label{eq22}
N_d(z) = \frac{N_d}{\sigma\sqrt{2\pi}} \exp\left[-\frac{(z - \mu)^2}{2\sigma^2}\right],$$

where

$$\label{eq23}
N_d = \frac{N N'}{\sum_{i = 1}^M N(z_i)}.$$

The correlation coefficient in Equation (\[eq17\]) is caculated using $N_d(z_i)$ instead of $N(z_i)$.

![\[fg02\] Example of discreet value of Gaussian distribution function generated by [gaussg]{} with $\mu = 0.5$ for $\sigma = 0.1$, $N' = 994$, $N_d = 109.665$ (solid line and square mark) and $\sigma = 0.04$, $N' = 998$, $N_d = 45.3348$ (dashed line and circle mark).](fig02.eps){width="9cm"}

![\[fg03\] Dependence of correlation coefficient $R^2$ on number of particles $N$ for $\mu = 0.5$, $\sigma = 0.1$, and $M = 10$.](fig03.eps){width="9cm"}

![\[fg04\] Dependence of correlation coefficient $R^2$ on number of groups $M$ for $\mu = 0.5$, $\sigma = 0.1$, and $N = 1000000$.](fig04.eps){width="9cm"}

Variation of number of particles $N$ and number of groups $M$ are also observed as illustrated in Figure \[fg03\] and Figure \[fg04\], respectively. It can be seen that larger $N$ gives better $R^2$ and larger $M$ gives bad $R^2$. Number of groups should be more than or equal to $N / M$ that the program [gaussg]{} can handled.

The next results are the sequences that produced from $N_d(z_i)$ as shown in Figure \[fg05\]. Only first four seeds are used to generate four sequences. These sequences has the same distribution function, which has $\mu = 0.5$,, $\sigma = 0.1$, $N = 100$, and $M = 10$. These results are produced by program [gausss]{}.

![\[fg05\] Sequences with seed: (a) 1, (b) 2, (c) 3, and (d) 4.](fig05a.eps "fig:"){width="12cm"}\
(a)\
![\[fg05\] Sequences with seed: (a) 1, (b) 2, (c) 3, and (d) 4.](fig05b.eps "fig:"){width="12cm"}\
(b)\
![\[fg05\] Sequences with seed: (a) 1, (b) 2, (c) 3, and (d) 4.](fig05c.eps "fig:"){width="12cm"}\
(c)\
![\[fg05\] Sequences with seed: (a) 1, (b) 2, (c) 3, and (d) 4.](fig05d.eps "fig:"){width="12cm"}\
(d)

Conclusion
==========

Two programs, [gaussg]{} for creating discreet groups and [gausss]{} for creating sequences, have been devoleped and tested. The discreen Gaussian distribution function can be produced. The sequences which has the same distribution function, can also be generated. Further investigation is needed how to register all available sequences for a distribution function. As $N$ increases the value $R^2$ approximates 1, but as $M$ increases the value $R^2$ decrease less than 1. $R^2 = 1$ can be achieved with larger $N$ and smaller $M$. The discreet Gaussian distribution function has different constant with its previously continuos distribution function which is used to generated the discreet and limited groups.

[**Acknowledgements**]{}\
\
Authors would like to thank to Institut Teknologi Bandung Alumni Association Research Grant in year 2011 for partially supporting to this work.

[**Appendix A: gaussg**]{}\
\

    /*
        gaussg.cpp
        Generate discreet groups of Gaussian distribution function
        Authors are Sparisoma Viridi and Veinardi Suendo
        Version date is 2011.07.17
    */

    #include <iostream>
    #include <fstream>
    #include <stdlib.h>
    #include <math.h>

    const double PI = 3.14159265;

    using namespace std;

    double Nz(double mu, double sigma, double N, double z);

    int main(int argc, char **argv) {
        if(argc < 6) {
            cout << "Version date is 2011.07.17" << endl;
            cout << "gaussg is written by Sparisoma Viridi ";
            cout << "and Veinardi Suendo" << endl;
            cout << "Generate discreet groups of particles ";
            cout << "that obey Gaussian distribution ";
            cout << "fuction" << endl;
            cout << endl;
            cout << "Usage: gaussg mu sigma N M output-file" << endl;
            cout << endl;
            cout << "All arguments are mandatory:" << endl;
            cout << "mu           average of Gaussian ";
            cout << "distribution function" << endl;
            cout << "sigma        width of Gaussian ";
            cout << "distribution function" << endl;
            cout << "N            number of particles" << endl;
            cout << "M            number of groups" << endl;
            cout << "output-file  output file" << endl;
        } else {
            double mu = atof(argv[1]);
            double sigma = atof(argv[2]);
            int N = atoi(argv[3]);
            int M = atoi(argv[4]);
            const char *ofn = argv[5];
            
            cout << "mu = " << mu << endl;
            cout << "sigma = " << sigma << endl;
            cout << "N = " << N << endl;
            cout << "M = " << M << endl;
            cout << "output-file = " << ofn << endl;
            
            double eps = 1E-3;
            double dz = mu * 1E-5;
            double NNz = N;
            
            double zmin = mu;
            while(NNz > eps) {
                NNz = Nz(mu, sigma, N, zmin);
                zmin -= dz;
            }
            cout << "zmin = " << zmin << endl;
            
            NNz = N;
            double zmax = mu;
            while(NNz > eps) {
                NNz = Nz(mu, sigma, N, zmax);
                zmax += dz;
            }
            cout << "zmax = " << zmax << endl;
            
            dz = (zmax - zmin) / M;
            cout << "dz = " << dz << endl;
            
            double zi[M];
            int Ni[M];
            double NN = 0;
            for(int i = 0; i < M; i++) {
                zi[i] = zmin + (i + 0.5) * dz;
                double z = zi[i];
                Ni[i] = (int) Nz(mu, sigma, N, z);
                NN += Ni[i];
            }
            cout << "N' = " << NN << endl;
            
            double NN2 = 0;
            for(int i = 0; i < M; i++) {
                Ni[i] = (int)(Ni[i] * (N/NN));
                // cout << i << "\t";
                // cout << Ni[i] << endl;
                NN2 += Ni[i];
            }
            cout << "N\" = " << NN2 << endl;
            
            double Nzi[M];
            
            double NN3 = 0;
            ofstream fout;
            fout.open(ofn);
            fout << "#i\tzi\tNi\tN(zi)" << endl;
            for(int i = 0; i < M; i++) {
                fout << i + 1 << "\t";
                fout << zi[i] << "\t";
                fout << Ni[i] << "\t";
                double z = zi[i];
                NN3 = 1.0 * N * NN2 / NN;
                Nzi[i] = Nz(mu, sigma, NN3, z);
                fout << Nzi[i] << endl;
            }
            fout.close();
            cout << "Nt = " << NN3 << endl;
            
            double SNi = 0;
            for(int i = 0; i < M; i++) {
                SNi += (Ni[i] * zi[i]);
            }
            double mui = SNi / NN2;
            
            double SStot = 0;
            double SSerr = 0;
            for(int i = 0; i < M; i++) {
                double dSStot = (Ni[i] - mui) * (Ni[i] - mui);
                SStot += dSStot;
                double dSSerr = (Ni[i] - Nzi[i]) * (Ni[i] - Nzi[i]);
                SSerr += dSSerr;
            }
            double R2 = 1 - SSerr/SStot;
            cout << "R^2 = " << R2 << endl;
        }
        return 0;
    }

    double Nz(double mu, double sigma, double N, double z) {
        double c1 = N  / (sigma * sqrt(2 * PI));
        double c2 = exp(-(z - mu)*(z - mu) / (2 * sigma * sigma));
        double c3 = c1 * c2;
        return c3;
    }

[**Appendix B: gausss**]{}\
\

    /*
        gausss.cpp
        Generate sequences from discreet groups of Gaussian
        distribution function
        Authors are Sparisoma Viridi and Veinardi Suendo
        Version date is 2011.07.17
    */

    #include <iostream>
    #include <fstream>
    #include <stdlib.h>
    #include <math.h>

    const double PI = 3.14159265;

    using namespace std;

    int main(int argc, char **argv) {
        if(argc < 4) {
            cout << "Version date is 2011.07.17" << endl;
            cout << "gausss is written by Sparisoma Viridi ";
            cout << "and Veinardi Suendo" << endl;
            cout << "Generate sequences from discreet groups ";
            cout << "of particles that\nobey Gaussian ";
            cout << "distribution fuction" << endl;
            cout << endl;
            cout << "Usage: seed input-file output-file" << endl;
            cout << endl;
            cout << "All arguments are mandatory:" << endl;
            cout << "seed         seed for random generator ";
            cout << "(1, 2, ..)";
            cout << endl;
            cout << "input-file   input file" << endl;
            cout << "output-file  output file" << endl;
        } else {
            long int seed = atoi(argv[1]);
            const char *ifn = argv[2];
            const char *ofn = argv[3];
            double mu = atof(argv[1]);
            
            cout << "seed = " << seed << endl;
            cout << "input-file = " << ifn << endl;
            cout << "output-file = " << ofn << endl;
            
            ifstream fin;
            fin.open(ifn);
            string buf;
            double d;
            long int i = 0;
            while(!fin.eof()) {
                fin >> buf;
                i++;
            }
            fin.close();
            
            int M = (int)((i - 1 - 4) / 4);
            double zi[M], Nzi[M];
            
            fin.open(ifn);
            int j = 0;
            fin >> buf; fin >> buf; fin >> buf; fin >> buf;
            while(!fin.eof() || j < i) {
                int k; fin >> k;
                fin >> zi[k-1];
                fin >> Nzi[k-1];
                fin >> buf;
                j++;
            }
            
            int N = 0;
            for(int l = 0; l < M; l++) {
                zi[l] = 0.001 * round(zi[l] * 1000);
                N += Nzi[l];
            }
            
            double seq0[N];
            int k = 0;
            for(int m = 0; m < M; m++) {
                for(int l = Nzi[m]; l > 0; l--) {
                    seq0[k] = zi[m];
                    k++;
                }
            }
            
            double seq1[N];
            for(int n = 0; n < N; n++) {
                seq1[n] = seq0[n];
            }
            
            srandom(seed);
            for(int n = 0; n < N; n++) {
                long int a1 = random();
                double a2 = 1.0 * a1 / RAND_MAX;
                int a3 = (int)(N * a2);
                swap(seq1[n], seq1[a3]);
            }
            
            ofstream fout;
            fout.open(ofn);
            for(int n = 0; n < N; n++) {
                fout << seq1[n] << endl;
            }
            fout.close();
        }
        return 0;
    }
