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 all
taillesim=5000;
for N=3:10
nb=0;
for k=1:taillesim
X=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);
end
tol=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;
end
end
probapprox(N-2)=nb/taillesim;
end
ns=[3:10];
for i=3:10
probhs(i-2)=probh(3,i);
end
scatter(ns,probhs)
hold on
plot(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
.
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 like
If 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 like
library(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(">=",N),"="),B)
So that I can solve here
and then, we simply have to check for the value of
slp$status
0 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, using
test=FALSE
x1=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(">=",N),"="),B)
if(slp$status==0) test=TRUE}
if(test==TRUE) nb0=nb0+1
So here, the code would be something like
taillesim=5000
probaap=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=0
for(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=FALSE
x1=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(">=",N),"="),B)
if(slp$status==0) test=TRUE
}
if(test==TRUE) nb0=nb0+1
}
probaap[N]=nb0/taillesim
probath[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!
Dans un 
(we focus only on solutions in [0,1]). Because here the game is symmetric, my son should also shoot at time 
, with density 
(we assume here that the density exists, or we seek only solution that are differentiable). Assume further that there exists 
>0 such that the support of my optimal shooting time (the time to shoot is now a random variable) is (
. Hence, if we differentiate (with respect to y), we have
.
, my expected payoff is
, i.e.
, where
,
, then 
=1, i.e. either 
=1 and then 
is constant, or 
is in proportional to 
.
=1/4 (but we don’t really care about that normalizing constant).
=0 if 
(duel started at time 0). Initially, the have no chance to kill the other one, and after one hour, they are certain to kill the other one. At time 
). What is the optimal time to shot ?
and 
for me (
stands for dad), and 
for the other player (here my son 
). The expected payoff becomes,
such that
, then 
and
. So
, then 
and
. So
. Hence, 
and 
are perfectly known, and the mixture parameter is the only one we care about.
(that cannot be observed), taking value when 
is drawn from 
and 
.
is known, denoted 
. Then I can predict the value of 
. And I can iterate from here.
is the best predictor of 
given my observations (as well as my belief in 
. Recall that we had 
was in 
, then we could have considered mean and standard deviations of observations such that 
before),
proportional to
, and 
being a parameter that will change, from 0 to 4.
or 
, depending whether 
or 
is the smallest parameter.
and 
is said to be an admissible solution for starting point x,
and 
. If we reformulate the dynamic programming idea, we obtain that if 
, then for all 
, sequence 
is a solution to problem 
. It comes that function v is a solution of Bellman’s equation
, 
and 
.
. The dynamic is then
. Assuming that
has a form that can be explicited).
.
.
.
est un ensemble convexe, compact. On parle alors de 
