Les basiques de Scilab

1. 👍 Soit $x$ un réel strictement positif et $n$ un entier naturel. On pose $g_n(x)=\dfrac{n!\times n^x}{x(x+1)...(x+n)}$. Écrire une script ou une fonction Scilab qui affiche ou retourne la valeurs de $g_n(x)$ pour $x>0$ et $n\in\N$ donnés. Écrire aussi un script qui trace la courbe représentative de $g_n$ sur $]0,1]$.

   👨‍🏫

   Voici une solution sous la forme d'une fonction:

   function r=g(n,x)
   	r=n^x
   	for k=0:n
   		r=r/(x+k)
   	end
   endfunction
   

   Pour la courbe de $g_n$:

   n=input('n=')
   x=linspace(0.1,1,100)
   for k=1:100
   	y(k)=g(n,x(k))
   end
   plot(x,y)
   

2. 👍 Montrer que, $\forall n\in\N^*$, l'équation "$ x\geq 0,\quad \dsum_{k=1}^n\dfrac{x^k}k=1$" possède une unique solution que l'on notera $x_n$, que la suite $(x_n)_{n\in\N^*}$ est strictement décroissante et possède une limite finie positive que l'on note $\ell$. Écrire un script Scilab qui calcule, par l'algorithme de dichotomie, une valeur approchée à $x_n$ à $\varepsilon$ près.

   👨‍🏫

   On vérifie que $f_n:x\mapsto \dsum_{k=1}^n\dfrac{x^k}k$ réalise une bijection strictement croissante de $\R^+$ sur lui même.

   On remarque que pour tout $n\in\N^*$, $x_n>0$ et $f_{n+1}(x_n)>f_n(x_n)$ i.e. $f_{n+1}(x_n) > 1$. D'où $x_{n+1} < x_n$.

   La suite $(x_n)_{n\in\N^*}$ est décroissante et minorée par $0$, elle est donc convergente.

   Le programme pour la fonction $f_ n$:

   function r=f(n,x)
   	r=0
   	for k=1:n
   		r=r+x^k/k
   	end
   endfunction
   

   On voit que $x_1=1$. La dichotomie:

   n=input('n=')
   vareps=input('epsilon=')
   a=0;b=1
   while (b-a)>vareps
   	c=(a+b)/2
   	if (f(n,a)-1)*(f(n,c)-1)<0
   		b=c
   	else
   		a=c
   	end
   end
   disp((a+b)/2)
   

3. 👍 Soit $\alpha$ un réel strictement positif. On considère une suite $(u_n)$ telle que $u_0>0$, $u_1>0$ et pour tout $n\geq 0$, $u_{n+2}=u_{n+1}+\alpha^n u_n$. On s'intéresse à la convergence de la suite $(u_n)$. Ecrire une fonction Scilab $\tt{ r=u(alpha,a,b,n)}$ qui renvoie le $n$-ième terme de la suite $u_n$ telle que $u_0=a$ et $u_1=b$. Utiliser cette fonction pour énoncer une conjecture sur la convergence de cette suite suivant la valeur de $\alpha$.

   👨‍🏫

   La fonction:

   function r=u(alpha,a,b,n)
   
   	if n<=1 
   		res=n*b+(1-n)*a
   	else
   		der=b;avtder=a
   		for k=2:n
   			tampon=der
   			der=der+avtder*alpha^(k-2)
   			avtder=tampon
   		end
   	end
   r=der
   endfunction
   

   On se doute que suivant que $\alpha$ est plus grand ou plus petit que $1$, le comportement asymptotique de $u_n$ risque d'être très différent.

   Le programme suivant permet de visualiser ce comportement (pour $\alpha=0.8,a=1,b=2$,$N=200$ et $\alpha=1.1,a=2,b=1$,$N=10$):

   
   alpha=input('alpha=')
   
   a=input('a=');b=input('b=')
   
   N=input('n max=')
   
   x=linespace(1,N,10)
   
   y=[]
   for k=1:10
   	y(k)=u(alpha,a,b,x(k))
   end
   
   plot(n,y,'*')
   

   On conjecture que la suite converge si $\alpha<1$ et diverge sinon.
4. 👍 Utiliser la méthode d'inversion pour réaliser une simulation avec Scilab d'une loi exponentielle de paramètre $1$. Utiliser cette fonction et la méthode de Monte-Carlo pour calculer une valeur approchée de $\Gamma(\alpha)$.

   👨‍🏫

   On sait que si $U$ suit la loi uniforme sur $]0,1]$ alors $X=-\ln(U)$ suit la loi $\Expo(1)$.

   On sait alors, d'après le théorème de transfert, que $X^{\alpha-1}$ admet une espérence qui vaut $\Gamma(\alpha)$. On peut alors appliquer la méthode de Monte-Carlo pour $\alpha>\frac 12$ car $X^{\alpha-1}$ admet alors une variance.

   Le programme:

   function r=expo(lambda)
   	r=-1/lambda*log(rand())
   end
   
   Esp=0
   alpha=input('alpha=')
   for k=1:100000
   	Esp=Esp+(expo(1))^(alpha-1)
   end
   disp(Esp/100000)
   

   Pour $\alpha=5$, on obtient environ $24.1$ pour $24$.
   Pour $\alpha\leq \frac 12$, on peut utiliser le fait que $\Gamma(\alpha+1)=\alpha\Gamma(\alpha)$.

5. 👍 On dit qu'une matrice carrée d'ordre $n$, $A=(a_{i,j})$, est à diagonale dominante si pour tout $i\in\zint{1,n}$, $|a_{i,i}|\geq \dsum_{j\neq i}|a_{i,j}|$. Écrire une fonction Scilab $\tt{estAdd(A)}$ qui retourne $\tt{\%T}$ si $A$ est à diagonale dominante et $\tt{\%F}$ sinon.

   👨‍🏫

   On remarque que la condition peut s'écrire, pour tout $i\in\zint{1,n}$, $2|a_{i,i}|\geq \dsum_{j=1}^n|a_{i,j}|$.

   D'où le programme:

   function rep=estAdd(A)
   	B=abs(A)
   	n=size(A)(1)
   	i=1
   	while i<=n & 2*B(i,i)>=sum(B(i,:))
   		i=i+1
   	end
   	rep=(i==(n+1))
   endfunction
   

6. 👍 On considère $X_{1,n},...,X_{r,n}$ des var qui suivent la loi uniforme discrète sur $\zint{1,n}$ et sont indépendantes. On pose $Y_{n,r}=\dfrac{\min(X_{1,n},...,X_{r,n})}n$ et $Z_{n,r}=\dfrac{\max(X_{1,n},...,X_{r,n})}{n}$. Ecrire deux fonctions Scilab $\tt{Min(n,r)}$ et $\tt{Max(n,r)}$ qui réalisent la simulation de $Y_{n,r}$ et de $Z_{n,r}$ et renvoient les valeurs simulées. En effectuant un grand nombre de simulations, représenter les densités empiriques de $Y_{100,r}$ et $Z_{100,r}$ sous la forme d'un histogramme comportant $20$ intervalles. Montrer la convergence en loi de ces deux suites et représenter les densités des lois limites. Étudier aussi la convergence en loi de $Y_{n,n}$ et $Z_{n,n}$.

   👨‍🏫

   • Simulations de $Y_{n,r}$ et $Z_{n,r}$:

     function res=Min(n,r)
     	res=min(grand(1,r,"uin",1,n))/n
     end
     function res=Max(n,r)
     	res=max(grand(1,r,"uin",1,n))/n
     end
     

   • Pour les histogrammes:

     r=input('r=')
     echantillonY=[]
     echantillonZ=[]
     
     for i=1:10000
     	echantillonY(i)=Min(100,r)
     	echantillonZ(i)=Max(100,r)
     end
     
     subplot(2,1,1)
     histplot((0:20)/20,echantillonY)
     subplot(2,1,2)
     histplot((0:20)/20,echantillonZ)
     

   • Cela donne:
     
   • Étudions la convergence en loi de $Y_{n,r}$. $Y_{n,r}$ est à valeurs dans $]0,1]$ et si $x\in]0,1]$ (on passe à ">" puisque c'est un min): $$\Pb([Y_{n,r}>x])=\Pb([\min(X_{1,n},...,X_{r,n})>nx])=\Pb(\bigcap_{1\leq k\leq r}[X_{k,n}>nx])$$ d'où par indépendance: $$\Pb([Y_{n,r}>x])\underset{\text{indép.}}=\prod_{k=1}^r\Pb([X_{k,n}>nx])=\prod_{k=1}^r(1-\Pb([X_{k,n}\leq nx])$$ donc finalement: $$\Pb([Y_{n,r}\leq x])=1-\left(1-\frac{\entiere{nx}}n\right)^r$$ et cette expression tend vers $1-(1-x)^r$ quand $n\to +\infty$. En résumé: $$\lim_{n\to +\infty}\Pb([Y_{n,r}\leq x])=\begin{cases} 0 & \text{ si }x\leq 0\\ 1-(1-x)^r & \text{ si }x\in]0,1]\\ 1 & \text{ si }x>1\end{cases}$$ On vérifie que la fonction limite est bien la fonction de répartition d'une variable aléatoire à densité. Une densité est par exemple la fonction $f_r$ définie par, pour tout réel $x$: $$f_r(x)=\begin{cases} 0 & \text{ si }x\notin [0,1]\\ r(1-x)^{r-1} & \text{ si }x\in[0,1]\end{cases}$$
   • Pour $Z_{n,r}$, on obtient: $$\lim_{n\to +\infty}\Pb([Z_{n,r}\leq x])=\begin{cases} 0 & \text{ si }x\leq 0\\ x^r & \text{ si }x\in]0,1]\\ 1 & \text{ si }x>1\end{cases}$$ On vérifie que la fonction limite est bien la fonction de répartition d'une variable aléatoire à densité. Une densité est par exemple la fonction $g_r$ définie par, pour tout réel $x$: $$g_r(x)=\begin{cases} 0 & \text{ si }x\notin [0,1]\\ rx^{r-1} & \text{ si }x\in[0,1]\end{cases}$$
   • Pour la représentation graphique:

     subplot(2,2,1)
     histplot((0:20)/20,echantillonY)
     subplot(2,2,2)
     histplot((0:20)/20,echantillonZ)
     x=linspace(0,1,100)
     y=r*(1-x).^(r-1)
     subplot(2,2,3)
     plot(x,y)
     z=r*x.^(r-1)
     subplot(2,2,4)
     plot(x,z)
     

     On obtient:
   • Si on s'intéresse à $Y_{n,n}$, on a pour $x\in]0,1[$: $$\Pb([Y_{n,n}\leq x])=1-\left(1-\frac{\entiere{nx}}n\right)^n=1-\e^{n\ln\left(1-\frac{\entiere{nx}}n\right)}$$ Or $\dlim_{n\to +\infty}n\ln\left(1-\frac{\entiere{nx}}n\right)=-\infty$, d'où pour tout $x\in]0,1[$: $$\dlim_{n\to +\infty}\Pb([Y_{n,n}\leq x])=1$$ De plus, si $x\leq 0$, $\Pb([Y_{n,n}\leq x])=0$ et lorsque $x\geq 1$, $\Pb([Y_{n,n}\leq x])=1$.

     On en conclut que $Y_{n,n}\cvl 0$ quand $n\to +\infty$.

     On montre de même que $Z_{n,n}\cvl 1$ quand $n\to +\infty$.

7. 👍 On considère une suite $(X_n)_{n\in\N^*}$ de variables indépendantes qui suivent la loi uniforme sur $\zint{1,N}$. Pour tout $n\in\N^*$, on note $Y_n$ le nombre de valeurs distinctes prises par $X_1,...,X_n$. Réaliser la simulation de $Y_n$ avec Scilab. Montrer que $Y_n$ est un estimateur asymptiotiquement sans biais et convergent de $N$.

   👨‍🏫

   • function res=Y(n,N)
     	res=0
     valeurs=1:N 
     for i=1:n
     	Xi=grand(1,1,"uin",1,N)
     	if valeurs(Xi)<>0
     		res=res+1
     		valeurs(Xi)=0
     	end
     end
     endfunction
     

   • $Y_n$ est une v.a.r. fonction qui ne dépend pas de $N$ de $(X_1,...,X_n)$. C'est donc bien un estimateur de $N$. Déterminons son espérance qui existe puisque $Y_n$ est bornée par $1$ et $N$.
   • On peut écrire que $Y_n=\dsum_{k=1}^N Z_k$, où $Z_k$ est la variable de Bernoulli égale à $1$ lorsqu'au moins l'un des $X_i$ prend la valeur $k$ et égale à $0$ sinon.
   • On a: $$\Pb([Z_k=0])=\Pb([X_1\neq k]\cap...\cap[X_n\neq k])\underset{\text{indép.}}=\Pb([X_1\neq k])\times...\times\Pb([X_n\neq k])=\left(1-\frac 1N\right)^n$$ D'où $\E(Z_k)=1-\left(1-\frac 1N\right)^n$.
   • Finalement, par linéarité, $\E(Y_n)=N\left(1-\left(1-\frac 1N\right)^n\right)$.
   • Quand $n\to +\infty$, $\left(1-\frac 1N\right)^n\to 0$, d'où $\E(Y_n)\to N$.
   • De plus on a $|Y_n-N|=N-Y_n$ d'où $\E(|Y_n-N|)=N-\E(Y_n)=N\left(1-\frac 1N\right)^n$. Ainsi par Markov, pour tout $\eps>0$: $$\Pb([|Y_n-N|\geq \eps])\leq \frac{N\left(1-\frac 1N\right)^n}{\eps}$$ et ce majorant tend vers $0$ quand $n\to +\infty$ ce qui prouve la convergence en probabilité de $Y_n$ vers $N$.
8. 😨 Soit $\B$ une base othonormée de $\R^n$, $F$ un sev de $\R^n$ dont $(u_1,...,u_p)$ est une base. On note $U=\M_{\B}(u_1,...,u_p)$, on note $U_{c,1}$,...,$U_{c,p}$ les colonnes de $U$ et pour tout $k\in\zint{1,p-1}$, $U_k$ la matrice formée des $k$ premières colonnes de $U$.
   Montrer que l'équation en $X_k$, $^tU_{k}U_{k}X_k=-^tU_{k}U_{c,k+1}$ admet une unique solution.

   💁‍♂️

   Montrer que $\rg(\,^tU_kU_k)=\rg(U_k)=k$ et en déduire qu'elle est inversible.

   
   Vérifier que la matrice $P$ dont les colonnes sont $U_{c,1},U_1X_1+U_{c,2},...,U_{p-1}X_{p-1}+U_{c,p}$ est la matrice dans la base $\B$ d'une base orthogonale de $F$.

   💁‍♂️

   Justifier assez rapidement que $\Vect(P_{c,1},...,P_{c,k})\subset\Vect(U_{c,1},...,U_{c,k})$ puis montrer l'égalité en prouvant par récurrence sur $k$ que $U_{c,k}\in\Vect(P_{c,1},...,P_{c,k})$.
   Pour finir, montrer que pour tout $k\in\zint{1,p-1}$, $\,^tU_kP_{c,k+1}=0$.

   Écrire une fonction Scilab $\tt{orthog(U)}$ qui renvoie la matrice $P$

   👨‍🏫

   • Pour tout $k\in\zint{1,p-1}$, $\rg(U_k)=k$ puisque $\rg(U)=p$ d'où comme nous l'avons vu à plusieurs reprises, $\rg(\,^tU_kU_k)=k$. Or cette matrice est carrée d'ordre $k$, elle est donc inversible.

     $-^tU_{k}C_{k+1}$ est une matrice colonne d'ordre $k$ d'où, il existe bien une unique matrice colonne d'ordre $k$, $X_k$ telle que $^tU_{k}U_{k}X_k=-^tU_{k}U_{c,k+1}$, c'est $(^tU_{k}U_{k})^{-1}(-^tU_{k}U_{c,k+1})$.

   • On note $P_{c,j}$ les colonnes de de $P$. On a clairement, pour tout $k\in\zint{1,p}$: $$\Vect(P_{c,1},...,P_{c,k})\subset\Vect(U_{c,1},...,U_{c,k})$$
   • On montre assez facilement par récurrence sur $k\in\zint{1,p}$ que: $$U_{c,k}\in\Vect(P_{c,1},...,P_{c,k})$$ Ceci prouve, en particulier, que $P$ est la matrice dans la base $\B$ d'une base de $F$.
   • Montrons que pour tout $k\in\zint{1,p-1}$, $P_{c,k+1}\in\Vect(U_{c,1},...,U_{c,k})^{\bot}$ ie $\,^tU_kP_{c,k+1}=0$.

     On a: $\,^tU_kP_{c,k+1}=\,^tU_k(U_kX_k+U_{c,k+1})=\,^tU_kU_kX_k+\,^tU_kU_{c,k+1}=0$ par définition de $X_k$.

   • Ainsi $P_{c,k+1}$ est orthogonal à $U_{c,1},...,U_{c,k}$ donc appartient à $\Vect(U_{c,1},...,U_{c,k})^{\bot}$ qui est aussi $\Vect(P_{c,1},...,P_{c,k})^{\bot}$, ce qui permet de conclure.
   • Le programme:

     function P=orthog2(U)
         p=size(U)(2) // nombre de colonne de U
         if rank(U)==p then
             P=zeros(U)
             P(:,1)=U(:,1) // colonne 1 de P = colonne 1 de U  
             for k=1:(p-1)
                 Uk=U(:,1:k) // matrice des k prem. colonnes de U 
                 Uc=U(:,k+1) // colonnes k+1 de U
                 P(:,k+1)=-Uk*(Uk'*Uk)^(-1)*Uk'*Uc+Uc
             end
         else
         P='la famille proposée n'est pas une base.'
         end
     endfunction
     

     Voici un exemple:

     --> P=orthog2([1 2 -1;2 1 0;0 3 0])
      P  = 
     
        1.   1.2  -0.6666667
        2.  -0.6   0.3333333
        0.   3.    0.3333333
     

     On vérifie que $\,^tPP$ est diagonale:

     --> P'*P
      ans  =
     
        5.         -2.220D-16   0.       
       -2.220D-16   10.8        0.       
        0.          0.          0.6666667
     

   .