{"word_count": 430, "character_count": 3719}
In my previous post, I tried to answer the following question Consider  points uniformly distributed on a sphere. What is the probability that the  points lie on a same hemisphere, for some hemisphere (there is no south or north here) ?If I have been able to use Monte Carlo simulations in dimension 2 (on a circle, not on a sphere), I could not get it in dimension 3. Hopefully, my colleague Simon gave me a nice solution (much more efficient than my previous code). Unfortunately, it was in Maple,clear alltaillesim=5000;for N=3:10nb=0; for k=1:taillesimX=randn(N,1);Y=randn(N,1);Z=randn(N,1);pts=[X Y Z];for i=1:N    pts(i,:)=pts(i,:)/norm(pts(i,:),2);endtol=1e-07;A=-pts;f=[0;0;0];b=zeros(N,1)-tol;%options=optimset('Display','off');[x,val,exitflag]=linprog(f,A,b);if(exitflag==1)  nb=nb+1;    endendprobapprox(N-2)=nb/taillesim;endns=[3:10];for i=3:10probhs(i-2)=probh(3,i);endscatter(ns,probhs)hold onplot(ns,probapprox);The idea is very clever (I did not know how to code it, but as we will see, it is actually simple - or say not too difficult - to code actually). The idea is based on the idea that all the points lie on the same hemisphere if there is some  such that  for all .&#160;That is not a big improvement, compared with the previous post, because we have to find such a vector. The idea suggested by Simon is to use some constraint optimization algorithm. We can try to optimize (maximize, or minimize, actually, we do not care) something likeIf the set of constraint is empty, then there is no solution, while we can find one (maybe more, but we do not care) if there is such a vector . I wanted to add some constraint on , like  (for the Euclidean norm), but that would be a non-linear constraint. So here, I chose to assume that  was such that . With R, the optimization routine would be likelibrary(lpSolve)A=rbind(pts,c(1,1,1))C=c(1,2,3)B=c(rep(1e-10,N),1)slp=lp("min",C,A,c(rep("&#62;=",N),"="),B)So that I can solve hereand then, we simply have to check for the value ofslp$status0 means that we found a vector  . But when you run the code, it does not work. Because actually, what is solved here is not exactly the program above, but the one below,(yes, it is mentioned in the help of the optimization function that only positive values are considered). So, a simple strategy, since we focus on an orthant here, is to consider all possible orthants. For instance, usingtest=FALSEx1=0:(2^3-1)%/%(2^(3-1))x2=((0:(2^3-1))-(2^2)*x1)%/%(2^(3-2))x3=((0:(2^3-1))-(2^2)*x1-2*x2)for(u in 1:(2^3)){pts2=cbind(pts[,1]*(-1)^x1[u],           pts[,2]*(-1)^x2[u],           pts[,3]*(-1)^x3[u])A=rbind(pts2,c(1,1,1))C=c(1,2,3)B=c(rep(1e-10,N),1)slp=lp("min",C,A,c(rep("&#62;=",N),"="),B)if(slp$status==0) test=TRUE}if(test==TRUE) nb0=nb0+1So here, the code would be something liketaillesim=5000probaap=rep(NA,10)probath=rep(NA,10)p=function(d,n) .5^(n-1) * sum(choose(n-1,0:(d-1)))for(N in 3:10){nb0=0for(k in 1:taillesim){MZ=matrix(rnorm(n*3),n,3)d=apply(MZ,1,function(z) sqrt(sum(z^2)))X=MZ[,1]/d; Y=MZ[,2]/d; Z=MZ[,3]/d; pts=cbind(X,Y,Z)test=FALSEx1=0:(2^3-1)%/%(2^(3-1))x2=((0:(2^3-1))-(2^2)*x1)%/%(2^(3-2))x3=((0:(2^3-1))-(2^2)*x1-2*x2)for(u in 1:(2^3)){pts2=cbind(pts[,1]*(-1)^x1[u],           pts[,2]*(-1)^x2[u],           pts[,3]*(-1)^x3[u])A=rbind(pts2,c(1,1,1))C=c(1,2,3)B=c(rep(1e-10,N),1)slp=lp("min",C,A,c(rep("&#62;=",N),"="),B)if(slp$status==0) test=TRUE}if(test==TRUE) nb0=nb0+1  }probaap[N]=nb0/taillesimprobath[N]=p(3,N)}and this time, the probability obtained using Monte Carlo simulation is extremely close to the theoretical value.And the algorithm is extremely fast! So using optimization routines to see if sets defined by linear constraints are empty - or not - is truly a great idea!