1 \chapter{Introduction}
3 As stated in \cite{scilabstats2001}, Scilab's core provides
4 a complete set of features related to simulation and statistical
5 computations.
6 Indeed, Scilab provides uniform pseudo-random number generators, functions to
7 compute the moments of a distribution and a complete set of
8 distributions. In this document, we will make a complete
9 overview on these features.
11 It must be noticed, though, that these features
12 are not as complete as in other languages, like R for example.
13 This is why several toolboxes have developped in order to extend
14 the features of Scilab. In this document, we will present two
15 major toolboxes, that is the Sci\_R toolbox and the Stix toolbox.
17 In the last chapter, we will analyse the missing statistical
18 features and will analyse how these features are available in other
19 tools, such as Matlab, R, or Octave.
21 \section{A sample session}
23 A good introduction on the statistical features of Scilab is \cite{scilabintro2007}.
24 In the remaining of this introduction chapter, we will try to have
25 a flavour of how to perform statistical computations with Scilab.
26 We focus on the algorithms which are used inside Scilab, to show what
27 exact algorithms perform the computations.
29 As a first example, we will generate a sequence of numbers from a
30 normal law with mean 0 and standard deviation 1 (example inspired and simplified
31 from \cite{scilabintro2007}). The probability density function (pdf)
32 and the cumulated probability density function of the normal law is
33 \begin{eqnarray}
34 f(x) &=& \frac{1}{\sqrt{2\pi}} e^{-\frac{t^2}{2}},\\
35 P(x) &=& \frac{1}{\sqrt{2\pi}} \int_{-\infty}^x e^{-\frac{t^2}{2}}.
36 \end{eqnarray}
38 The empirical cumulated density function \cite{artcomputerKnuthVol2} of a given
39 set of data $\{x_i\}_{i=1,N}$ is given by
40 \begin{eqnarray}
41 F_N(x) &=& \frac{\textrm{number of } x_1,x_2,\ldots,x_n \textrm{ that are }\leq x}{N}.
42 \end{eqnarray}
45 The numerical method used by Scilab to generate such numbers is the Polar
46 method for normal deviates, as presented in \cite{artcomputerKnuthVol2}.
48 \lstset{language=Scilab}
49 \lstset{numbers=left}
50 \lstset{basicstyle=\footnotesize}
51 \lstset{keywordstyle=\color{green}\bfseries}
52 \begin{lstlisting}
53 N=200;
54 x = rand(1,N,"normal");
55 Xsorted =gsort(x,"g","i");
56 Ydata = (1:N)/N;
57 plot(Xsorted,Ydata);
58 e=gce();
59 e.children.polyline_style=2;
60 xtitle("Empirical Cumulated Density Function of Normal Law with 200 samples")
61 filename = "introduction_ecdfnormal.png";
62 xs2png(0,filename);
63 \end{lstlisting}
65 The empirical cumulated density function
66 is presented in figure \ref{introduction_ecdfnormal}.
68 \begin{figure}[htbp]
69 \begin{center}
70 \includegraphics[height=10cm]{introduction_ecdfnormal.png}
71 \end{center}
72 \caption{Empirical Cumulated Density Function of Normal Law with 200 samples}
73 \label{introduction_ecdfnormal}
74 \end{figure}
76 To compare the data which is produced by rand with the
77 cumulated density function of the normal law, we use the
78 \emph{cdfnor} primitive. This primitive is based on \cite{Algorithm715}
79 and uses rational functions that theoretically approximate the normal
80 distribution function to at least 18 significant decimal digits. The same
81 primitive can be used to compute the inverse of the cumulated density
82 function. In that case, rational functions are used as starting values to
83 Newton's Iterations which compute the inverse standard normal.
85 \lstset{language=Scilab}
86 \lstset{numbers=left}
87 \lstset{basicstyle=\footnotesize}
88 \lstset{keywordstyle=\color{green}\bfseries}
89 \begin{lstlisting}
90 N=200;
91 x = rand(1,N,"normal");
92 Xsorted =gsort(x,"g","i");
93 Ydata = (1:N)/N;
94 x=linspace(-3,3,100);
95 P=cdfnor("PQ",x,zeros(x),ones(x));
96 plot(Xsorted,Ydata,x,P);
97 \end{lstlisting}
99 The comparison plot between the empirical cdf and the
100 computed cdf is presented in figure \ref{introduction_ecdcomparison}.
102 \begin{figure}[htbp]
103 \begin{center}
104 \includegraphics[height=10cm]{introduction_ecdfcomparison.png}
105 \end{center}
106 \caption{Cumulated Density Function of Normal Law : comparison of cdf from
107 rational functions and empirical cdf from Polar method }
108 \label{introduction_ecdcomparison}
109 \end{figure}
111 The moments of a distribution can be computed with the
112 \emph{mean}, \emph{variance} and \emph{stdev} Scilab macros,
113 which are implementations of the moments. For the variance
114 and standard deviation, the scaling factor is $N-1$.
115 In the following script, one computes these moment for
116 an increasing number of samples, from $10^1$ to $10^5$.
118 \lstset{language=Scilab}
119 \lstset{numbers=left}
120 \lstset{basicstyle=\footnotesize}
121 \lstset{keywordstyle=\color{green}\bfseries}
122 \begin{lstlisting}
123 nbpoints = 5;
124 means=zeros(nbpoints,1);
125 vars=zeros(nbpoints,1);
126 stdevs=zeros(nbpoints,1);
127 nlist = 1:nbpoints;
128 for i = nlist
129   N=10^i;
130   x = rand(1,N,"normal");
131   means(i) = mean(x);
132   vars(i) = variance(x);
133   stdevs(i) = stdev(x);
134 end
135 plot(nlist,[means,vars,stdevs]);
136 \end{lstlisting}
138 The convergence plot of the moments is presented in
139 figure \ref{introduction_convergencemoments}.
141 \begin{figure}[htbp]
142 \begin{center}
143 \includegraphics[height=10cm]{introduction_convergencemoments.png}
144 \end{center}
145 \caption{Convergence of the moments of the normal law}
146 \label{introduction_convergencemoments}
147 \end{figure}