Tag Archives: géométrie

Enveloppe convexe de points tirés au hasard

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 ?

Être – en même temps – moyen et extraordinaire

Juste avant 1850, Adolphe Quételet publiait Sur l’appréciation des documents statistiques, et en particulier sur l’application des moyennes, et introdusait alors le concept d’homme moyen (l’ouvrage est cité dans La théorie des moyennes et son emploi dans les sciences d’observation, de Joseph Lottin, paru en 1909 – je n’ai pas trouvé en ligne cet ouvrage de Quételet). Ce concept sera détaillé encore davantage dans Anthropométrie, ou mesure des différentes facultés de l’homme, concept qui sera abondament illustré

Mais si la moyenne est un concept fondamental en statistique, il convient de garder en mémoire que la moyenne n’existe peut être pas, au sein de la population. C’est ce que rappelle la fameuse ‘blague’ “un mec qui à la tête dans le four et les pieds dans le frigo, en moyenne, il est bien“. La moyenne peut être quelquechose d’extraordinaire (au sens premier du terme).

Continue reading Être – en même temps – moyen et extraordinaire

Pavage et magie, rotations et symétries

(oui, mon fils voulait une rime dans le titre… promis on va tenter de faire une comédie musicale à partir de mes billets). La semaine passée, un article intéressant était publié (http://www.bbc.co.uk/news/…) sur sciences et magie. Cet article m’a fait penser à une réflexion que je m’étais faite lors de notre visite au MoMaths, cet hiver. Pour revenir sur l’histoire que j’avais raconté il y a quelques mois, sur la fin de la visite, les enfants avaient joué avec les pavages, à la Escher (les pièces étaient magnétiques, et s’assemblaient sur un grand tableau blanc)

Le soucis c’est qu’à aucun moment je n’ai pu voir d’explication du  pourquoi ça marche ? ou pour être plus explicite comment a été construit le singe pour que tout s’emboîte aussi bien. Car il y a un truc… La magie ça n’existe pas ! ce sont juste des maths (c’est un peu ce que racontent Persi Diaconis et Ron Graham dans Magical Mathematics, même si j’avoue ne pas l’avoir fini, faute de temps, et d’avoir laissé des chapitres en attente). Mon fils travaillant les rotations et les symétries, à l’école, j’ai voulu qu’on prenne le temps de revenir sur les pavages. Mais avant de me lancer avec lui, j’ai essayé de m’entraîner un peu tout seul (dimanche après midi, pendant le temps calme).

La difficulté (ou la magie) est qu’il y a des symétries dans tous les sens. Regardons un exemple simple (ou disons pas trop compliqué) avec des chinois (identifiés par les chapeaux dits chinois). En regardant rapidement, on note que le triangle ci-dessous est reproduit à l’infini, par un jeu de rotations. Il suffit donc d’avoir un tel triangle. Le soucis est que le découpage (parfaitement symétrique) n’est pas trivial (il est délicat de faire un personnage en faisant un découpage dans un triangle, même si c’est la méthode retenue dans l’activité mentionnée sur https://irem.univ-lille1.fr/…)

Le secret part de l’observation suivante: sur le graphique de gauche, en dessous, les trois points (qui sont des intersections des trois personnages, ou disons des trois couleurs) sont les sommets d’un triangle équilatéral. C’est là que la magie va opérer, car c’est cette figure (avec une autre, juxtaposée) que l’on va tenter de reproduire. En fait, un personnage (disons un petit chinois bleu) c’est deux triangles collés (on doit pouvoir parler de parallélogramme isocèle).

On ne le voit pas très bien, mais pourtant, au départ, on avait un parallélogramme isocèle (dont les contours sont tracés en noir) qui était  tout bleu. Bien entendu, on a découpé des parties et en les a envoyé ailleurs (par rotation autour d’un sommet). Le plus simple est de construire la figure, à la main. On part donc du parallélogramme suivant, tout bleu :

Oui, j’ai repris les fonctions de base utilisées l’autre jour pour expliquer R aux enfants,

source("http://freakonometrics.free.fr/RforKIDZ.R")

En fait, j’ai souhaité définir une fonction créant un polynôme à partir d’une liste de points (ou de sommets).

polyg=function(L,couleur="blanc",tour=couleur,ep=1,st=1,...){
  b=unlist(L)
  n=length(b)
  polygon(b[seq(1,n,by=2)],b[seq(2,n,by=2)],col=couleurc(couleur),border=couleurc(tour),lwd=ep,lty=st,...)}

La fonction pour faire une rotation d’une liste de points est ici

rotation=function(L,angle=90,centre=L[[1]]){
  b=unlist(L)
  n=length(b)
  X=cbind(b[seq(1,n,by=2)]-centre[1],b[seq(2,n,by=2)]-centre[2])
  Y=matrix(c(cos(angle/180*pi),sin(angle/180*pi),-sin(angle/180*pi),cos(angle/180*pi)),2,2)%*%t(X)
  Z=Y+centre
  L2=L
  for(i in 1:(n/2)){L2[[i]]=Z[,i]}
  return(L2)}

On a alors ici les quatre sommets

A=c(4,0)
B=c(4,10)
C=c(4+sqrt(3)*5,5)
D=c(4+sqrt(3)*5,15)

et on peut tracer le parallélogramme avec

L1=list(A,B,D,C)
polyg(L1,"bleu",tour="noir",ep=3)

On peut alors jouer à découper des parties, puis à les pivoter (on va faire des rotations, quoi). Par exemple, le prend le triangle en haut à droite, et je pivote de 240 degrés autour d’un des sommets du parallélogramme. C’est la figure ci-dessus à droite.  Le triangle, il est défini par le code suivant

E=c(10,10)
F=c(7,10+sqrt(3)) 
L1=list(B,E,F)
polyg(L1,tour="rouge",ep=3)

ce qui enlève la portion – en traçant un triangle blanc – et pour faire la rotation, j’utilise le code suivant

L2=rotation(L1,240)
polyg(L2,"bleu","rouge",ep=3)

(la seconde ligne permet de tracer le triangle bleu). Il va de soi qu’on peut découper des objets plus complexes que des triangles, voire non connexes. Par exemple, je peux découper la forme verte ci-dessous, puis la pivoter par rapport au sommet en haut à gauche, de 240 degrés.

G=c(10,2*sqrt(3)+10)
H=c(4+sqrt(3)*5,2*sqrt(3)+10)
F1=c(8,10+2/sqrt(3))
F2=c(8,10+4/sqrt(3))
L1=list(B,E,F1,F2,B,G,H,D,G)
polyg(L1,"vert",tour=NA)
L2=rotation(L1,240)
polygL(L2,"bleu",tour=NA)

(il n’est pas nécessaire de revenir au point B, mais peu importe). Ou la figure jaune, qui va pivoter autour du sommet en bas à droite, de 120 degrés.

F3=c(8,10-2/sqrt(3))
F4=c(8,10-4/sqrt(3))
G2=c(10,-2*sqrt(3)+10)
H2=c(4+sqrt(3)*5,-2*sqrt(3)+10)
L1=list(H,G,F2,F1,E,F3,F4,G2,H2,H)
polyg(L1,"jaune",tour=NA)
L2=rotation(L1,120,centre=C)
polygL(L2,"bleu",tour=NA)

Avec des ciseaux, ça sera plus simple, mais le but est aussi de voir que programmer, c’est facile. En gros, on a deux fonctions de base: une qui trace un polynôme et l’autre qui fait des rotations.

Si on fait les deux en même temps (à peu près), on enlève toute la partie en haut à droite, et on la répartie: une partie à gauche, et une partie en bas à droite.  On a alors une forme qui ressemble à un chinois. Non ? On peut aller creuser un peu sous les bras, pour améliorer les mains, mais on y est presque, me semble-t-il…

On peut continuer à améliorer, par exemple en enlevant des parties, puis en faisant des rotations par rapport à un des sommets initiaux. Comme le petit triangle jaune, qui va pivoter de 120 degrés par rapport au sommet en bas à droite. Puis le triangle vert par rapport au sommet en bas à gauche (on note qu’on pivote par rapport aux sommets du triangle équilatéral qu’on avait repéré, initialement)

S=c(8-.4,-2*sqrt(3)+10-.4)
L1=list(F4,G2,S,F4)
polyg(L1,"jaune")
L2=rotation(L1,120,centre=C)
polyg(L2,"bleu")

Bon, là on fait une pause, parce qu’on a, par morceaux, notre chinois ! On un paquet de lettre (voire un paquet de sommets qu’il va falloir nommer, suite à nos rotations). On les nomme avec un code qui ressemble à (oui, on aurait pu pu les nommer dès le départ, lors des rotations)

rE=rotation(list(E),240,centre=B)[[1]]
rF1=rotation(list(F1),240,centre=B)[[1]]
rH=rotation(list(H),240,centre=B)[[1]]

Une fois qu’on a identifié nos sommets, on a notre chinois,

chinois=list(B,rE,rF1,LG,S2,LF2,rH,A,r2H,r2G,r2S2,
        r2F2,r2F1,r2E,rF3,rF4,r2S,rG2,rH2,C,H2,G2,S,F4,F3,E,B)

Bon, ensuite, reste à paver, ce qui est est le plus simple. Par exemple, on peut faire des rotations,

r1chinois=rotation(chinois,120,centre=A)
polyg(r1chinois,couleur="jaune")

voire définir une translation,

translation=function(L,h){
b=unlist(L)
n=length(b)
X=cbind(b[seq(1,n,by=2)]+h[1],b[seq(2,n,by=2)]++h[2])
L2=L
for(i in 1:(n/2)){L2[[i]]=X[i,]}
return(L2)}

puis translater notre chinois,

r1chinois=rotation(chinois,120,centre=A)
polyg(r1chinois,couleur="jaune")
polyg(translation(r1chinois,15/sqrt(3),15),couleur="jaune")

Il manque l’autre, le vert. On peut vérifier qu’il se superpose parfaitement :

Bon, ben ça y est, on a constitué notre pavage. Reste à travailler un peu le code, essentiellement sur la forme, et on pourra jouer avec les enfants, à faire les lézards, des singes, et toutes sortes de pavages possibles, à la Escher. En notant qu’on peut utiliser d’autres formes de bases plus complexes, comme expliqué sur http://xavier.hubaut.info/…

Régression médiane et géométrie

Suite à mon exposé d’hier, Pierre me demandait d’argumenter un point que j’avais évoqué oralement sans le justifier: “ une régression médiane passe forcément par deux points du nuage “. Par exemple,

x=c(1,2,3)
y=c(3,7,8)
plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
library(quantreg)
abline(rq(y~x,tau=.5),col="red")

Essayons de le justifier… de manière un peu heuristique. Commençons un peu au hasard, avec une droite qui passe “entre” les points (on admettra que la régression médiane sépare l’espace en deux, avec la moitié des points en dessous, la moitié au dessus, à un près pour des histoires de parité).

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=-1,b=3.2,col="blue")

c’est joli, mais on peut sûrement faire mieux. En particulier, on va chercher ici à minimiser la somme des valeurs absolue des erreurs (c’est le principe de la régression médiane). Essayons de translater la courbe, vers le haut, ou vers le bas,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=-1,b=3.2,col="blue",lwd=2)
for(i in seq(-2,3,by=.5)){
abline(a=i,b=3.2,col="blue",lty=2)}

Si on regarde ce que vaut la somme des valeurs absolues des erreurs, on a

d=function(h) sum(abs(y-(h+3.2*x)))
H=seq(-4,6,by=.01)
D=Vectorize(d)(H)
plot(H,D,type="l")

Formellement, l’optimum est ici

> optimize(d,lower=-5,upper=5)
$minimum
[1] -0.2000091

$objective
[1] 2.200009

Bref, retenons cette courbe, et notons qu’elle passe par un des points,

Et c’est assez normal. On commençait avec deux points au dessus, et un en dessous. Soit la somme des valeurs absolues des erreurs . Si on translate de , on passe de à (tant que l’on a toujours deux points au dessus, car celui en dessous sera toujours en dessous). Si on translate de , (là encore tant qu’il n’y a qu’un point en dessous). Bref, on a intérêt à translater vers le haut. Si on dépasse le premier point rencontré, on se retrouve dans la situation inverse, avec un point au dessus et deux au dessous, et on a intérêt à redescendre. Bref, la translation optimale revient a s’arrêter dès qu’on croise un point. Autrement dit, la régression passe forcément par un point. Au moins.

Maintenant, essayons de faire pivoter la courbe autour de ce point, là encore afin de minimiser la somme des valeurs absolues des erreurs,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=optimize(d,lower=-5,upper=5)$minimum,b=3.2,
col="blue",lwd=2)
points(x[1],y[1],cex=1.8,lwd=2,col="blue")
for(i in seq(-1,5,by=.25)){
abline(a=(y[1]-i*x[1]),b=i,col="blue",lty=2)}

La distance est alors, en fonction de la pente de la droite

d2=function(h) sum(abs(y-((y[1]-h*x[1])+h*x)))
H=seq(-4,6,by=.01)
D=Vectorize(d2)(H)
plot(H,D,type="l")

Là encore on peut formaliser un peu

> optimize(d2,lower=-5,upper=5)
$minimum
[1] 2.500018

$objective
[1] 1.500018

Et si on regarde cette dernière droite, on passe par deux points,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
h=optimize(d2,lower=-5,upper=5)$minimum
abline(a=(y[1]-h*x[1]),
b=h,col="purple",lwd=2)

Pour comprendre pourquoi, comparons les deux cas,

  • en faisant un pivot vers le bas, pour un des points, la valeur absolue des erreurs augmente alors que pour l’autre, elle diminue

  • en faisant un pivot vers le haut, là encore pour un des points, la valeur absolue des erreurs augmente alors que pour l’autre, elle diminue, mais en sens inverse par rapport au cas précédent,

Et pour faire simple, dans le premier cas, le gain (sur la baisse) compense la perte, alors que dans le second cas, c’est le contraire. Bref, on a intérêt à pivoter vers le bas, ou vers le point à droite. Jusqu’à l’atteindre. Et optimalement, on passera par ce point. Bref, on passera par deux points. Et je laisse les plus courageux regarder avec plus de trois points, mais c’est toujours pareil…

Discussions géométriques du soir

J’avais raconté l’autre jour l’aventure de mon fils qui souhaitait faire des triangles avec des bouts de bois (ici). Il a récidivé ce soir. Sauf qu’il n’avait qu’un grand baton, et que j’avais fermé le portillon qui l’empêchait d’aller dans la forêt en chercher (la nuance entre romancer et devenir mythomane doit être mince). Je lui ai suggéré de le casser en trois portions, puis de constituer un triangle avec les trois morceaux restant. Manque de bol, l’expérience rate (comme l’autre fois). Et me souvenant de mes pénibles calculs de la dernière fois, je me suis rappelé qu’il avait une chance sur deux de pouvoir effectivement construire un triangle… sauf que la manière d’obtenir les trois bouts de bois me semble différente… donc peut être que ce 1/2 n’est plus valide ?
En fait, je crois qu’il faut définir clairement comment on va couper les bouts de bois.

  • méthode 1: le coupage simultané

Je suppose que je peux couper d’un coup le bout de bois en deux. Formellement, je tire deux points sur l’intevalle correspondant à mon bout de bois, les tirages étant uniformes, et indépendants. On arrive alors à trois morceaux, alatex2png.2.php.pngalatex2png-1.2.php.png et alatex2png-2.2.php.pngen sachant que 

alatex2png-3.2.php.png

la longueur du bout de bois dont je diposais initialement. On supposera que aalatex2png.2.png=1 pour simplifier (ce qui ne changera rien au résultat final, on s’en doute). L’ensemble des points que je peux obtenir peuvent alors se visualiser comme des points du simplexe, 

alatex2png-4.2.php.png

autrement dit, on peut les représenter dans le triangle équilatéral ci-dessous.
Or on ne peut constituer un triangle que si la plus grande des longeurs n’excède pas la somme des deux autres. Supposons que alatex2png-2.2.php.png soit la plus grande longueur (et qu’elle soit visualisée graphiquement par rapport à la distance au côté en bas), le point rouge ci-dessous (correpondant à un triplet de longueurs possibles) ne convient pas. En fait, dès lors que l’on entre dans l’aire rosée, on ne peut plus constituer un triangle…

Moralité, à moins de tomber dans le triangle jaune, on ne peut pas constituer de triangle…. Ce qui veut dire qu’on n’a plus qu’une chance sur quatre !

  • méthode 2: le coupage séquentiel

Une autre méthodde consister à couper d’abord en deux (en tirant un point au hasard sur mon intervalle), puis de recouper en deux le plus grand morceau… Damned, mon fils a de bonne idée, mais ça sent les calculs à n’en plus finir ça…Ou alors va falloir continuer à essayer de faire des dessins, car mine de rien, c’est plus simple que les calculs de la dernière fois. Bon, la longueur de mon bout de bois est la hauteur de mon triangle. Je commence par le couper en deux, la distance au côté en bas devant être la plus petite des deux.

Une fois fixé la taille de ce petit bout, il me reste à couper le plus grand morceaux en deux. Pour cela, je vais me déplacer le long du segment rouge ci-dessous,
La position sur ce petit segment me donnant ensuite les deux autres longueurs des morceaux. Or rappelons que pour que nos bouts de bois forments un triangle, il faut forcément appartenir au triangle jaune,
Pour résumer, nous avions exclus la partie supérieure du triangle, et conditionnellement au fait que nous sommes dans le lozange inférieur, nous souhaitons être dans le triangle jaune. Arg, cette fois-ci la probabilité passe à une chance sur trois !

  • couper d’un coup comme découpage séquentiel aléatoire ?

Mais attends deux minutes mon bonhomme (oui, je parle à mon fils). Si j’avais choisi de découper le plus petit morceau, on sait qu’il est impossible de constituer un triangle car par définition alatex2png-5.2.php.png Supposons que l’on casse le baton en deux, et que l’on tire ensuite au hasard (avec probabilité 1/2 et 1/2) quelle partie sera alors coupée en deux, on doit se ramener au cas précédant non ? La probabilité de constituer un triangle doit alors être

alatex2png-6.2.php.png

Damned, mais ça ne fait pas un quart ça ! Bon, je crois que je vais prendre une aspirine et aller demander à mon fils d’aller ranger ses bâtons !

Discussions géométriques du matin (partie 1)

Petit résumé du débat du matin avec mon fils …
– papa, papa j’ai trouvé trois bouts de bois,
– ben… fais une figure géométrique avec *
– ça marche pas !
– comment ça, “ça marche pas” ?!?
– j’arrive pas à faire un triangle !
– … hein ? attends, j’arrive
genre “mais qu’est-ce que tu racontes ? t’es pas réveillé ce matin ou quoi ?“. Bref, après avoir vu qu’effectivement ça n’était pas possible (un des bouts de bois était beaucoup plus long que les deux autres),  je me suis demandé si mon fils n’avait vraiment pas eu de chance… ou si au contraire, c’était plutôt “normal” de ne pas pouvoir faire un triangle à partir de trois bouts de bois….
Soient https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-01.pnghttps://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-02.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-03.png les longueurs des trois morceaux de bois (oui, une fois le café avalé, je peux prendre un peu de temps pour formaliser). Ces trois bouts de bois forment  un triangle si et seulement si la somme des deux plus petits est plus grande que le plus grand. Autrement dit, si je m’intéresse aux statistiques d’ordre

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-04.png

il faut que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-05.png

Pour faire simple, on supposera que les https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-50.png suivent une loi uniforme sur [0,1].

  • Approche par simulation

Par simulation, ça ne devrait pas être trop compliqué d’avoir un ordre de grandeur de la probabilité

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-06.png

Le code est tout simplement

> n=1000000
> U=matrix(runif(3*n),n,3)
> U1=apply(U,1,max)
> U2=apply(U,1,sum)-U1
> U3=U1<U2
> mean(U3)
[1] 0.501067

Finalement la probabilité semble être proche de 1/2… ça doit pouvoir se montrer….

  • Approche par des calculs de proba

Il doit être possible de calculer directement

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-06.png

ou bien, on peut aussi noter que

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-07.png

ce qui ne fait intervenir que deux variables (certes, toujours non indépendantes, mais on devrait s’en sortir malgré tout): le maximum des trois variables, et la somme des trois variables. On cherche alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-08.png

Je pense qu’il est plus simple de travailler sur la première expression, en notant qu’il existe 6 permutations de ces 3 valeurs, et cette probabilité s’écrit simplement

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-10.png

(qui  est un des 6 arrangements possibles, tous équiprobables). Cette probabilité est alors tout simplement, compte tenu de l’indépendance des tailles des bouts de bois, ce qui est peut être une hypothèse un peu forte, certes,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-11.png

qui doit pouvoir se réécrire

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-12.png

L’intégrale au centre se réécrit alors simplement, à condition de savoir comment https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-13.png se positionne par rapport à 1. Deux cas assez différents sont à distinguer: comme on veut https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-134.png, soit https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-15.png, soit https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-16.png. Autrement dit on écrira l’ensemble sur lequel on calcule la première intégrale comme la réunion de deux ensembles (disjoints),

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-18.png

et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-19.png

Sur le premier ensemble,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-20.png

et sur le second

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-21.pngBref, l’intégrale se retrouve sous la forme de deux intégrales

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-22.png

et

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-23.png

La première intégrale correspond à un cas où https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-24.png est petit, ou plutôt https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-25.png. La seconde intégrale est un cas où https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-26.png est grand, ou plutôt https://perso.univ-rennes1.fr/arthur.charpentier/latex/b<br />
ois-27.png. Bref, la première intégrale peut s’écrire

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-28.png

et la seconde

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-29.png

(on aura en mémoire le facteur 6, que je ne mes pas pour alléger un peu). Bref, on a deux fois la même intégrale, qui s’écrire (dans le second cas à une translation des indices près)

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-30.png

Notre probabilité est alors 6 fois le double de cette intégrale, soit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/bois-31.png

Bon, finalement c’est beaucoup plus improbable de pouvoir construire un triangle que je ne le croyais au départ…. Damned ! Moralité, mieux vaut ne pas trop réfléchir avant d’avoir pris un café…. Et comme “papa” est très fort en mauvaise foi, “papa” expliquera que si la longueur des bois de bois n’était pas distribuée suivant une loi uniforme, l’intuition initiale n’était pas complètement stupide… Non mais !
* pour ma défense, c’était ça ou la révision de la table de multiplication de 5… en gros, “papa a pas encore pris son café, c’est pas le moment“.