Typo
[scilab.git] / scilab_doc / statisticswithscilab / introduction.tex
1 \chapter{Introduction}
2
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.
10
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.
16
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.
20
21 \section{A sample session}
22
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.
28
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}
37
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}
43
44
45 The numerical method used by Scilab to generate such numbers is the Polar 
46 method for normal deviates, as presented in \cite{artcomputerKnuthVol2}.
47
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}
64
65 The empirical cumulated density function 
66 is presented in figure \ref{introduction_ecdfnormal}.
67
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}
75
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.
84
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}
98
99 The comparison plot between the empirical cdf and the 
100 computed cdf is presented in figure \ref{introduction_ecdcomparison}.
101
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}
110
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$.
117
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}
137
138 The convergence plot of the moments is presented in 
139 figure \ref{introduction_convergencemoments}.
140
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}
148
149