Volume of n-sphere? Please help me with my version:
4 次查看(过去 30 天)
显示 更早的评论
Hello Experts,
I have to write a function for 10-dimensional sphere using Monte - Carlo method.
What I did is this:
N = 1000000;
g = 1;
S = 0;
RandVars = -1 + 2*rand(1,10,N);
for i = 1:N
if (sum(RandVars(1,10,:).^2) < 1)
S = S+1;
end
end
Please revise my code and assist me with computing the volume using the Monte Carlo.
Many thanks, Steve
0 个评论
采纳的回答
John D'Errico
2013-11-11
My guess is this homework was chosen to show you how inefficient Monte Carlo can be when used in a high number of dimensions. THe curse of dimensionality strikes here.
The known volume of a unit 10-sphere is simple enough to compute. Google is your friend here, but I'll just use a tool I've posted on the FEX.
format long g
spheresegmentvolume([],10)
ans =
2.55016403987735
Use of Monte Carlo to compute volume is a dartboard approach. Thus throw a series darts at a hypercube of edge length 2 that just contains the unit 10-sphere. Count the number of times the dart falls inside the sphere, and divide by the total number of darts thrown. Multiply by the known 10-cube volume, and you get an approximation of the 10-sphere volume.
The trick is to do this efficiently. Do NOT use a loop.
% the dimension of the problem
dim = 10
% the volume of a 10-cube of edge length 2
cubevol = 2^dim;
% sample size
N = 10e6;
% draw the samples
X = 2*rand(N,dim) - 1;
% Distance from the origin for each sample
% Note that I could have skipped the square root, since this is a unit sphere.
R = sqrt(sum(X.^2,2));
% count the number of hits
C = sum(R <= 1)
C =
24925
% compute the estimated volume
V = C/N*cubevol
V =
2.55232
As you can see, the computed volume was not terribly inaccurate compared to the known volume of 2.5502… 3 significant digits is as much as we can expect.
To be honest, if I were your instructor, I'd have specified a 20 or even 50 dimensional sphere, to make that hit rate vanishingly small. For the 10-sphere, the fraction of darts that actually hit the sphere is only about 0.25%. But that is a large enough fraction that we got almost 3 digits of accuracy from our pretty large sample size.
hitrate = spheresegmentvolume([],10)/cubevol
hitrate =
0.00249039457019272
A nice thing about Monte Carlo is is allows you to compute a measure of your uncertainty in the final estimate.
A binomial random variable with p = 0.00249039… has variance n*p*(1-p). So if we wish to put +/- 2*sigma confidence limits around the number of hits we should have expected:
C + [-2 2]*sqrt(N*hitrate*(1-hitrate))
ans =
24609.7735731207 25240.2264268793
And of course, the final volume will have confidence bounds of:
V + cubevol/N*[-2 2]*sqrt(N*hitrate*(1-hitrate))
ans =
2.52004081388756 2.58459918611244
Again, learn how to avoid loops when you can do so.
6 个评论
Marc
2013-11-12
Does this mean John is out of retirement?? Or are you simply taking pity on poor Steve? Either way it is good to see your name in this user forum.
John D'Errico
2013-11-12
@Marc - I can never truly retire because there are constantly things I need/want to maintain on the tools I've posted here. A few neat things are on the horizon too. (Major VPI upgrade! SLM gui tool one day...) So I look in occasionally, and if I see something on which I can add something meaningful, I might chime in.
更多回答(1 个)
Walter Roberson
2013-11-10
To calculate volume, you need to divide the successes, S, by the number of attempts, to get the fill fraction. Then multiply the fill fraction by the area of a hypercube each side of which is the same as the "2" you used in "2 * rand(1,100,N)"
Note, the leading 1 in your rand() is not serving any purpose.
In each iteration of your loop, you are summing the same thing. Where are you indexing by "i" ?
另请参阅
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!