Results 1 to 3 of 3

Thread: Monte Carlo calculation of volumes of n-balls

  1. #1
    Member dcromley's Avatar
    Join Date
    Apr 2009
    Location
    Wyoming, USA
    Age
    86
    Posts
    80
    Rep Power
    23

    Monte Carlo calculation of volumes of n-balls

    I could have sworn that somebody started a thread here regarding an article in The American Scientist.
    HTML Code:
    http://www.americanscientist.org/issues/pub/an-adventure-in-the-nth-dimension/1
    So I looked at it, got interested, and wanted to verify some of the "zaniness" of the subject. He is working with "n-balls" of radius 1, which are "spheres" in n dimensions. A 3-ball is the normal sphere in 3 dimensions; the volume is 4/3pi*r^2 = 4.1888. A "2-ball" is a circle in 2 dimensions; the "volume" is 2*pi*r = 3.1416. A 1-ball is a line from -1 to 1; the "volume" is 2. Going the other way -- going up in dimensions is not clear at all, but he gives a table on page 4 and a formula on page 8.

    My thought was to write a program to divide the "n-cube" into 0.01 increments and count the points that were "inside" the "n-ball" by comparing the distance from the center to 1. Oops -- in 10 dimensions, that would be 100^10 = 100,000,000,000,000,000,000 points. Maybe increments of 0.1. That's still 10^10 = 10,000,000,000 points.

    Another approach is the "Monte Carlo" method
    HTML Code:
    http://en.wikipedia.org/wiki/Monte_Carlo_method
    which is picking many random points and counting the ones that are inside the "sphere". The distance from the center is the square root of the sum of the squares of the n coordinates. The volume is the ratio of points "inside" to total points times the volume of the "n-cube".

    This TB program does that, and surprise, surprise -- the results agree with the table! The higher dimension "n-balls" are slow in converging because the "volume" of the ball is very small compared to the vast "volume" of the "n-cube" of size 2. For dimension 10, the "n-ball" is 2.55 and the "n-cube" is 1024 cubic units.

    Changing the program to size 1 (n-ball in a "unit" n-cube") gives more normal results -- the volumes don't increase.

    Table of volumes of n-balls
    Dim  Size 2  Size 1
    1    2.0000  1.0000 
    2    3.1416  0.7854 
    3    4.1888  0.5236 
    4    4.9348  0.3084 
    5    5.2638  0.1645 
    6    5.1677  0.0807 
    7    4.7248  0.0369 
    8    4.0587  0.0159 
    9    3.2985  0.0064 
    10   2.5502  0.0025
    
      Uses "Console", "UI"
      Global hWnd As Long
      hWnd = Canvas_Window "    " & _ 
        "Monte Carlo Calculation of Volumes of ""N-Balls"" in size 2 ""N-Cubes"" in Dimensions 1-10;" & _ 
        "    ESC To Exit", 0, 0, 1018, 736
      ' x goes from 0-1017; y from 0-735
      ' dCromley, Dec, 2011
      Const ymax As Number = 6 ' for y scaling
      Const maxdim As Number = 10 ' max dimensions
      Canvas_Attach hWnd, 0, %TRUE
      Canvas_Scale(0,ymax, 1017,0)
      Canvas_Clear Rgb(255,255,255)
      Global aIns(maxdim), aIters(maxdim) As Long, aVol(maxdim) As Single
      Global canvasx, niterations As Long ' ix is screenx
    
    Sub TBMain()
    ' http://www.americanscientist.org/issues/pub/an-adventure-in-the-nth-dimension/
      Local i1, i2, idim As Long, rsqr, vol As Single
      Randomize
      ygrid()
      For canvasx = 1 To 1017
        niterations = canvasx^2 ' 1 to 1000000 or so
        For idim = 1 To maxdim ' for each of 10 dimensions
          vol = GetVol(idim) ' volume is ratio of ins / ins+outs
          vol *= 2^idim ' adjust for size 2 sides
          aVol(idim) = vol ' save for printing
          Canvas_SetPixel(canvasx, vol)      
        Next idim                   
        If canvasx > 100 And Mod(canvasx, 10) = 0 Then zPrint() ' every once in awhile
        If GetAsyncKeyState(%VK_ESCAPE) Then Exit Sub
      Next canvasx         
      Beep
      Canvas_WaitKey
    End Sub                         
    
    Function GetVol(nDim As Long) As Single
    ' -- get random points in a n-dimension box 
    '    return vol, the fraction of points within ball
      Local Ins, iIter, iDim As Long, rsqr, vol As Single
      For iIter = aIters(ndim)+1 To nIterations
        Incr aIters(ndim) 
        rsqr = 0
        For iDim = 1 To nDim ' add up the iDim sides^2
          rsqr += (Rnd-.5)^2 ' -.5 to +.5 squared
        Next iDim           
        ' don't take sqrroot, just comare with .25
        If rSqr < .25 Then Incr aIns(ndim) ' if within ball, incr Ins
      Next iIter                             
      vol = aIns(ndim)/aIters(ndim) ' volume is fraction of Ins / Ins+Outs
      Function = vol
    End Function
    
    Sub zPrint() ' labels and values
      Local idim, ix As Long
      For idim = 1 To maxdim     
        ix = canvasx - idim*10
        If Mod(ix, 100) = 0 Then 
          Canvas_SetPos(canvasx, aVol(idim))
          Canvas_Print idim ' label
        End If  
        Canvas_SetPos(80+idim*80, .8)
        Canvas_Print idim & ": " & Format$(aVol(idim), "0.0000")
      Next idim                     
      Canvas_SetPos(50, .8)
      Canvas_Print "Iters="&niterations ' # iterations
      Canvas_Redraw
    End Sub
    
    Sub ygrid()
      Local iy As Long
      For iy = 0 To ymax-1
        Canvas_SetPos(0, iy+.1)
        Canvas_Print Format$(iy, "0.00")
        Canvas_Line((0, iy), (1017, iy), Rgb(192,192,192))
      Next iy
    End Sub
    
    Last edited by dcromley; 11-12-2011 at 17:59. Reason: 3 to 4 digits after decimal

  2. #2
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    Here is a similar problem which I solved with thinBasic, using the Monte Carlo method.

    Dan

    http://www.thinbasic.com/community/s...8473#post78473

    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  3. #3
    Member dcromley's Avatar
    Join Date
    Apr 2009
    Location
    Wyoming, USA
    Age
    86
    Posts
    80
    Rep Power
    23
    Very interesting also, thanks.

    As I was running mine, I knew how dependent it was on "rnd". And reading your and Eros' dialog on random number generation was interesting. Quite a bit of analysis to that. I'll look for "Numerical Recipes in C". Regards, Dave

Similar Threads

  1. calculation failure?
    By largo_winch in forum thinBasic General
    Replies: 1
    Last Post: 19-07-2011, 13:46

Members who have read this thread: 0

There are no members to list at the moment.

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •