Le week-end dernier, Jean-Baptiste qui était de passage à la maison, me présentait un problème amusant de géométrie, lié à un papier mis en ligne l’an dernier, monotonicity of facet numbers of random convex hulls. Dans cet article, ils montrent que quand on tire n points au hasard (dans un espace de dimension d) alors P_n, le nombre moyen de face de l’enveloppe convexe est strictement croissant avec n. Si on n’a pas regardé la démonstration (on avait mieux à faire), Jean-Baptiste me disait que ce problème très simple était en fait très complexe. Et bien entendu, comme ça m’a interpelé, j’ai voulu regarder plus en détails. En dimension d=2, et en tirant en plus uniformément sur le carré unité (oui, j’ai fait très très simple).

Pour tirer des points au hasard, et récupérer l’enveloppe convexe, c’est assez simple, par exemple avec 15 points

library(sp)
library(geosphere)
n=15
UV=matrix(runif(2*n),n,2)
CH=chull(UV)
PLCH=UV[c(CH,CH[1]),]
plot(c(0,0,1,1),c(1,0,0,1))
polygon(PLCH,border="blue",col=rgb(0,0,1,.3))
points(UV,pch=19,cex=2,col="red")

On le voit, ici l’enveloppe convexe est un polygône de 8 côtés (ou 8 sommets). On peut alors faire un petit code pour tirer des points au hasard, construire l’enveloppe convexe, et sortir des infos (nombre de points extrémaux, surface de l’enveloppe convexe, présence ou non de certains points – sur la diagonale, etc)

simu=function(n,isplot=FALSE){
UV=matrix(runif(2*n),n,2)
CH=chull(UV)
PLCH=UV[c(CH,CH[1]),]
nb_ex=length(CH)
p_in=function(u) point.in.polygon(u,u,PLCH[,1],PLCH[,2])
pts_in=Vectorize(p_in)(seq(.5,.95,by=.05))
if(isplot==TRUE) lines(PLCH,col=rgb(0,0,1,.25))
return(list(nb=nb_ex,area=areaPolygon(PLCH),pts=pts_in))}

par exemple

plot(c(0,0,1,1),c(1,0,0,1),col="white")
for(s in 1:1000){S=simu(5,isplot=TRUE)}

(on pourrait bien entendu stocker tout plein de choses)

ou encore avec n= 20 points au lieu de 5

plot(c(0,0,1,1),c(1,0,0,1),col="white")
for(s in 1:1000){S=simu(20,isplot=TRUE)}

Essayons de boucler un peu sur n maintenant

Np=c(3,4,5,6,7,8,10,15,20,30,40,50,75,100,200)
VN=VA=rep(NA,15)
NN=matrix(NA,20000,15)
VPT=matrix(NA,15,10)
for(i in 1:15){
N=A=rep(NA,20000)
PT=matrix(NA,20000,10)
np=Np[i]
for(s in 1:20000){
S=simu(np,isplot=FALSE)
N[s]=S$nb
PT[s,]=S$pts
A[s]=S$area
}
NN[,i]=N
VN[i]=mean(N,na.rm=TRUE)
VA[i]=mean(A,na.rm=TRUE)
VPT[i,]=apply(PT,2,function(x) mean(x,na.rm=TRUE))
}

Cette fois on stocke tout plein de choses. On peut juste faire un boxplot du nombre des points extrémaux en fonction de la taille de l’échantillon

VV=rep(Np,each=20000)
boxplot(as.vector(NN)~as.factor(VV))

Oui, en moyenne, ça semble croitre. Plus amusant, si on regarde la moyenne en fonction de \log(n)

plot(Np,VN,type="l",log="x",col="blue")

on obtient… une belle droite ! Le nombre moyen de points extrémaux croit en \log(n). On peut même avoir la pente de cette droite,

> lm(VN~log(Np))
Coefficients:
(Intercept) log(Np)
0.05224 2.58717717

Je laisse les plus courageux trouver du sens à ce 2.58717… Si on continue un peu, on peut regarder la probabilité que (u,u) soit à l’intérieur, pour plusieurs valeurs de u.

plot(Np,VPT[,10],type="l",log="x",col="blue")
lines(Np,VPT[,8])
lines(Np,VPT[,5],col="red",lwd=2)

On retrouve des fonctions croisantes, en fonction de n, mais la convexité semble dépendre de l’endroit où se trouve le point u. Amusant, non ?