\section{Simulations}

\subsection{R package}

The developed R-package is called \verb"morpheus" \cite{Loum_Auder} and divided into two main parts:
\begin{enumerate}
	\item the computation of the directions matrix $\mu$, based on the empirical
		cross-moments as described in the previous sections;
	\item the optimization of all parameters (including $\mu$), using the initially estimated
		directions as a starting point.
\end{enumerate}
The former is a straightforward translation of the mathematical formulas (file R/computeMu.R),
while the latter calls R constrOptim() method on the objective function expression and its
derivative (file R/optimParams.R). For usage examples, please refer to the package help.

\subsection{Experiments}
In this section, we evaluate our algorithm in a first step using mean squared error (MSE). In a second step, we compare experimentally our moments method (morpheus package \cite{Loum_Auder}) and the likelihood method (with felxmix package \cite{bg-papers:Gruen+Leisch:2007a}). 
We arbitrarily choose the parameters for the simulations, which should be discovered by the algorithms (ours, and the likelihood algorithm).

Experiment 1 (dimension 2):
\begin{align*}
	K &= 2\\
	p &= (0.5, 0.5)\\
	b &= (-0.2, 0.5)\\
	\beta &=
		\begin{pmatrix}
		1 & 3\\
		-2 & 1
		\end{pmatrix}
\end{align*}

Experiment 2 (dimension 5):
\begin{align*}
	K &= 2\\
	p &= (0.5, 0.5)\\
	b &= (-0.2, 0.5)\\
	\beta &=
		\begin{pmatrix}
			1 & 2\\
			2 & -3\\
			-1 & 0\\
			0 & 1\\
			3 & 0
		\end{pmatrix}
\end{align*}

Experiment 3 (dimension 10):
\begin{align*}
	K &= 3\\
	p &= (0.3, 0.3, 0.4)\\
	b &= (-0.2, 0, 0.5)\\
	\beta &=
		\begin{pmatrix}
			1 & 2 & -1\\
			2 & -3 & 1\\
			-1 & 0 & 3\\
			0 & 1 & -1\\
			3 & 0 & 0\\
			4 & -1 & 0\\
			-1 & -4 & 2\\
			-3 & 3 & 0\\
			0 & 2 & 1\\
			2 & 0 & -2
		\end{pmatrix}
\end{align*}

For all three experiments we use both logit and probit links.
Computations are always run on the same data both for our package and flexmix $-$ which is a reference for this kind of estimation, using an iterative algorithm to maximize the log-likelihood. Results are aggregated over $N=1000$ Monte-Carlo runs.\\


\noindent \textbf{Mean squared error (MSE)}.
Graphical representations of the MSE versus the sample size ($n$) are given in figures~\ref{mse_2}, ~\ref{mse_5} and \ref{mse_10}. In each figure, we represent the MSE associated with each parameter vector versus the sample size. We can see that the goodness of the estimation depends on  the sample size and that enough observations is needed to properly estimate the parameters. We have from figure~\ref{mse_2} and figure~\ref{mse_5} that for our moment method, with dimension less or equal to $5,$ the necessary sample size is around of $10^5.$ For large dimension, figures~\ref{mse_10} show that $10^6$ is not enough to estimate the parameters vectors $\beta.$ Indeed our estimation method use the spectral estimator, which need more enough data, as starting point.\\






\begin{figure}[!!ht!!]
\centering
	\includegraphics[keepaspectratio=true,height=3.5cm,width=8cm]{logit_2.png}\\
	\includegraphics[keepaspectratio=true,height=3.5cm,width=8cm]{probit_2.png}
	\caption{Experiment 1; Top: \textit{logit} link function $\Big($(a) MSE($\widehat{\beta}$), (b) MSE($\widehat{b}$,$\widehat{p}$)$\Big)$, bottom: \textit{probit} link function: $\Big($(c) MSE$(\widehat{\beta})$, (d) MSE($(\widehat{b}$,$\widehat{p}$)$\Big)$.}
	\label{mse_2}
\end{figure}

\begin{figure}[!!ht!!]
\centering
	\includegraphics[keepaspectratio=true,height=3.5cm,width=8cm]{logit_5.png}\\
	\includegraphics[keepaspectratio=true,height=3.5cm,width=8cm]{probit_5.png}
	\caption{Experiment 2; Top: \textit{logit} link function $\Big($(a) MSE($\widehat{\beta}$), (b) MSE($\widehat{b}$,$\widehat{p}$)$\Big)$, bottom: \textit{probit} link function: $\Big($ (c) MSE$(\widehat{\beta})$, (d) MSE($(\widehat{b}$,$\widehat{p}$)$\Big)$.}
	\label{mse_5}
\end{figure}


\begin{figure}[!!ht!!]
\centering
	\includegraphics[keepaspectratio=true,height=3.5cm,width=8cm]{logit_10.png}\\
	\includegraphics[keepaspectratio=true,height=3.5cm,width=8cm]{probit_10.png}
	\caption{Experiment 3; Top: \textit{logit} link function $\Big($(a) MSE($\widehat{\beta}$), (b) MSE($\widehat{b}$,$\widehat{p}$)$\Big)$, bottom: \textit{probit} link function: $\Big($ (c) MSE$(\widehat{\beta})$, (d) MSE($(\widehat{b}$,$\widehat{p}$)$\Big)$.}
	\label{mse_10}
\end{figure}

%\begin{figure}[!!ht!!]
%\centering
%	\includegraphics[keepaspectratio=true,height=4cm,width=8cm]{Figures/logit_2.png}
%	\caption{Experiment 1: \textit{logit} link function for our algorithm; (a) MSE($\widehat{\beta}$), (b) MSE($\widehat{b}$,$\widehat{p}$).}
%	\label{mse_logit_2}
%\end{figure}
%
%\begin{figure}[!!ht!!]
%\centering
%	\includegraphics[keepaspectratio=true,height=4cm,width=8cm]{Figures/logit_5.png}
%	\caption{Experiment 2: \textit{logit} link functionfor our algorithm; (a) MSE($\widehat{\beta}$), (b) MSE($\widehat{b}$,$\widehat{p}$).}
%	\label{mse_logit_5}
%\end{figure}
%
%
%\begin{figure}[!!ht!!]
%\centering
%	\includegraphics[keepaspectratio=true,height=4cm,width=8cm]{Figures/logit_10.png}
%	\caption{Experiment 3: \textit{logit} link function for our algorithm; (a) MSE($\widehat{\beta}$), (b) MSE($\widehat{b}$,$\widehat{p}$).}
%	\label{mse_logit_10}
%\end{figure}
%
%\begin{figure}[!!ht!!]
%\centering
%	\includegraphics[keepaspectratio=true,height=4cm,width=8cm]{Figures/probit_2.png}
%	\caption{Experiment 1: \textit{probit} link function for our algorithm; (a) MSE($\widehat{\beta}$), (b) MSE($\widehat{b}$,$\widehat{p}$).}
%	\label{mse_probit_2}
%\end{figure}
%
%\begin{figure}[!!ht!!]
%\centering
%	\includegraphics[keepaspectratio=true,height=4cm,width=8cm]{Figures/probit_5.png}
%	\caption{Experiment 2: \textit{probit} link function for our algorithm; (a) MSE($\widehat{\beta}$), (b) MSE($\widehat{b}$,$\widehat{p}$).}
%	\label{mse_probit_5}
%\end{figure}
%
%\begin{figure}[!!ht!!]
%\centering
%	\includegraphics[keepaspectratio=true,height=4cm,width=8cm]{Figures/probit_10.png}
%	\caption{Experiment 3: \textit{probit} link function for our algorithm; (a) MSE($\widehat{\beta}$), (b) MSE($\widehat{b}$,$\widehat{p}$).}
%	\label{mse_probit_10}
%\end{figure}





\noindent \textbf{Algorithms performance}.
To evaluate algorithms performance, the total number of sample points is fixed to $n = 10^5.$ This value still enough to observe correct performances, yet small enough to remain realistic.

On the figures~\ref{res_logit} and \ref{res_probit}, all (true) parameters ($p, b, \beta$) are re-ordered in a real vector, of size $K \times (d+2) - 1$. This vector is plotted as a line to improve visualization experience, but it must be noted that this does not represent any curve data. Dotted lines corresponds to the computed values plus or minus one standard deviation, and computed values themselves are represented with long dashed lines. The leftmost column corresponds to experiment 1 ($d=2$), the middle one to experiment 2 ($d=5$), and the rightmost column corresponds to experiment 3 ($d=10$).




\noindent \textit{Figure \ref{res_logit}: logit link.}
While most of the times flexmix finds a better solution than our proposed algorithm (smaller variance), both methods are good on average for $d \leq 5$.
The case $d=10$ is not handled well neither by the flexmix package nor by our package (the latter showing even poorer accuracy).
Indeed in this relatively high dimension the number of observations should be much higher. 

\begin{figure}[!!ht!!]
\centering
	\includegraphics[keepaspectratio=true,height=8cm,width=12cm]{res_logit.png}
	\caption{\textit{Logit} link function. Top: our package, bottom: flexmix. From left to right: experiment 1, 2 and 3 respectively}
	\label{res_logit}
\end{figure}
\noindent \textit{Figure \ref{res_probit}: probit link.}
In this case our algorithm performs slightly better than its flexmix counterpart for $d \leq 5$.
However, again, the variance in the case $d=10$ is way too high $-$ in fact even the average value is generally wrong, when coefficients are non-zero.
We can increase $n$ by a factor $100$ to obtain more accurate results.\\
\begin{figure}[!!ht!!]
\centering
	\includegraphics[keepaspectratio=true,height=8cm,width=12cm]{res_probit.png}
	\caption{\textit{Probit} link function. Top: our package, bottom: flexmix. From left to right: experiment 1, 2 and 3 respectively}
	\label{res_probit}
\end{figure}



% TODO: improve placement: see https://tex.stackexchange.com/questions/8625/force-figure-placement-in-text

Considering both links with $d \leq 5$, the algorithms performances are comparable with a global small advantage to our method.
The case $d=10$ is handled better by the flexmix algorithm, although clearly not well.
Finally, concerning our method we observe a tendancy to overestimate large parameters while underestimating small ones.
This observation is inversed for flexmix.\\

\noindent \textbf{Computational time}.
It must be noted that our algorithm timing does not depend on $n$, since it operates on matrices of size $O(d \times K)$.
Thus it can be more suitable for very large datasets, where the variability is clearly reduced. Figure~\ref{timing} (time by seconds versus $\log_{10}n$) illustrate this fact: the timings clearly favor our package $-$ because increasing $n$
has almost no impact on the running time. However, flexmix timings are not that high:
just a few minutes for the longest run, on average, on one million sample points.
To obtain the data shown on the figure we averaged $1000$ runs on random parameters.\\

\begin{figure}
\centering
	\includegraphics[height=8cm,width=10cm]{timings.png}
	\caption{Timings versus sample size. Dotted lines: flexmix; others: our algorithm. Black for $d=2,$ red for $d=5$ and blue for $d=10.$}
	\label{timing}
\end{figure}


%\begin{figure}[h!]
%\centering
%	\includegraphics[height=6cm,width=14cm]{Figures/timings_logit.png}
%	\caption{\textit{Logit} link function. Timings versus sample size. Our algorithm on the left, flexmix on the right. $d=2$ in long-dashed line, $d=5$ in dotted line, $d=10$ in solid line.}
%	\label{timing_logit}
%\end{figure}
%
%
%\begin{figure}[h!]
%\centering
%	\includegraphics[height=6cm,width=14cm]{Figures/timings_probit.png}
%	\caption{\textit{Probit} link function. Timings versus sample size. Our algorithm on the left, flexmix on the right. $d=2$ in long-dashed line, $d=5$ in dotted line, $d=10$ in solid line.}
%	\label{timing_probit}
%\end{figure}