[Insight-users] random number generator

[Stefan Klein]

Hi Jim,

You observed well. I think this is the expected behaviour. I do not think 
the truncation issues are the cause.

M = the number of pixels in the image
N = the number of loops over the image (200 in the example)
R = The number of pixels selected by the RandomIterator (= 100M in this 
example)

image2:
after one loop through the image:
         P( pixel_i = 0x sampled ) = 1/2      for i = 1...M
         P( pixel_i = 1x sampled ) = 1/2 for i = 1...M
after N loops through the image:
         P (pixel_i = K times sampled) = \sum_{k=1...N} binom_coeff(N,K) 
(1/2)^K (1/2)^{N-K}     for i = 1...M
         (which is of course the binomial distribution)
         Expectation: E( K ) = N/2 = 200/2 = 100
         Variance: Var( K ) = N/4 = 200/4 = 50
         Standard Deviation: SD( K ) = sqrt (N/4) = 7. as you observed.

image 4:
after one jump of the RandomIterator
         P( pixel_i = 0x sampled ) = (M-1)/M       for i = 1...M
         P( pixel_i = 1x sampled ) = 1/M        for i = 1...M
after R jumps of the RandomIterator:
         P( pixel_i = K times sampled ) =  \sum_{k=1...R} binom_coeff(R,K) 
(1/M)^K ((M-1)/M)^{R-K}      for i = 1...M
         (so again a binomial distribution)
         Expectation: E(K) = R/M = 100M / M = 100
         Variance: Var(K) = R/M * (M-1)/M = approximately R/M (if M is big) 
= 100.
         Standard Deviation: SD(K) = sqrt( var(K) ) = 10. as you observed.

My statistic math classes were some time ago, so correct me if i made a 
mistake :)

Regards,
Stefan.

At 16:43 16/03/05, Miller, James V (Research) wrote:
>An interesting observation:
>
>The standard deviation of the pixels in image2 and image4 are suprisingly 
>different.  image2 has a standard deviation of 7 and image4 has a standard 
>deviation of 10.
>
>If I understand the original post correctly, both images were created 
>using the same random number generator.  image2 was created by drawing a 
>random number 200 times at each pixel and counting the number of times the 
>random number was greater than 0.5.  image4 was created using a 
>ImageRandomIteratorWithIndex to visit pixels in the image at random, 
>performing enough random visits so that each pixel was on average visited 
>100 times (the same expected value as the pixels in image2).
>
>Naively, I would expect the mean and standard deviations of these two 
>images to be very close.
>
>The RandomIterator pretty much just draws a random number between 0 and 
>the number of pixels and unrolls the resulting linear index into an index 
>in the image.  The truncation of the random number to an integer and then 
>to an index increases the standard deviation of how many times a pixel is 
>visited from 7 to 10.
>
>Probably not a useful observation.... but interesting.
>
>Jim
>-----Original Message-----
>From: insight-developers-bounces at itk.org 
>[mailto:insight-developers-bounces at itk.org]On Behalf Of Blezek, Daniel J 
>(Research)
>Sent: Wednesday, March 16, 2005 10:05 AM
>To: Stefan Klein; Insight-users at itk.org; insight-developers at itk.org
>Subject: [Insight-developers] RE: [Insight-users] random number generator
>
>Stefan,  I started a thread on this a while ago relating to 
>non-deterministic behavior across platforms.  If the ITK design committee 
>approves the move, it would be great to standardize the random number 
>generator to be sufficiently random, fast and generate the same sequences 
>(from a given seed) across platforms.  The discussion culminated with 
>Brad's post: 
><http://www.itk.org/mailman/private/insight-developers/2005-January/006220.html>http://www.itk.org/mailman/private/insight-developers/2005-January/006220.html 
>I'm not sure any action happened on the suggestions.
>
>As you point out, many registration algorithms depend on random sampling, 
>exactly where I came across the problem.
>
>Cheers,
>-dan
>-----Original Message-----
>From: insight-users-bounces at itk.org 
>[mailto:insight-users-bounces at itk.org]On Behalf Of Stefan Klein
>Sent: Wednesday, March 16, 2005 9:16 AM
>To: Insight-users at itk.org
>Subject: [Insight-users] random number generator
>
>
>Dear itk-users,
>
>I did some tests with the underlying random generator of the 
>itkImageRandomIteratorWithIndex and it seems that, in Windows, it is 'not 
>very random'.
>
>The itkImageRandomIterator uses the following random-generator:
><ITKSOURCE>\Utilities\vxl\core\vnl\vnl_sample.<h/cxx>
>
>In Linux the drand48 random-generator is used, which is a good choice.
>In Windows however a "simple congruential random number generator" is 
>implemented, since drand48 is not available in Windows. This gives 
>inferior results, in my experience.
>
>To get a feel for the result look at the image1.<mhd/raw>, which is in the 
>file: randomtestresults.zip, which you can download from: 
>http://www.isi.uu.nl/People/Stefan/
>The gray-values in this image show how many times a pixel was sampled. The 
>sampling
>process works was defined as follows:
>    "An ImageRegionIterator walks N=200 times through the image and tests
>    at each voxel whether to sample it or not. The test is performed by
>    drawing a number between 0 and 1 using the random generator defined in
>    vnl_sample.h; a value >=0.5 means that the voxel is sampled"
>As you can see in the image the pixels are not really selected at random...
>
>If we use the itkImageRandomIteratorWithIndex to sample the image, the 
>result (image3.<mhd,raw>) may look better at first sight, but the 
>histogram looks terrible.
>
>On the internet I found the following link:
>http://paine.wiau.man.ac.uk/pub/doc_vxl/core/vnl/html/classvnl__random.html
>It describes the files vnl_random.<h,cxx>; those files are not included in 
>the vxl-version that comes with ITK. I tried to use these as a random 
>generator. The results are now much better. Image2, which is generated in 
>the same was as image1, but with the random generator defined in 
>vnl_random, does not show any structure anymore. Image4, which is 
>generated using an ImageRandomIterator that uses the vnl_random as its 
>underlying random number generator, has the histogram that you would expect.
>
>For more details please look at the source of the test program, which you 
>can also download from http://www.isi.uu.nl/People/Stefan/ : 
>itkrandomtestsource.zip. You may reproduce the results with this program. 
>Note that only in Windows bad results will be obtained. In Linux there is 
>no problem.
>
>Algorithms that rely on a good random generator may fail if you use the 
>itkRandomImageIteratorWithIndex under Windows. For example, in my research 
>on nonrigid registration I tried to use a stochastic gradient descent 
>optimisation method for minimising the MattesMutualInformation in 
>combination with a B-spline transform (it may speed up your registration 
>algorithms; if you are interested:
>http://www.isi.uu.nl/Research/Publications/publicationview.php?id=1011 ). 
>When I used the itkImageRandomConstIteratorWithIndex for selecting spatial 
>samples the registration results got significantly worse than when I used 
>the itkImageMoreRandomConstIteratorWithIndex (which uses the vnl_random as 
>underlying random number generator).
>
>Any comments on this would be appreciated!
>
>Stefan
>
>
>
>
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20050317/c71bb98c/attachment-0001.htm