[Insight-users] random number generator

[Miller, James V (Research)]

Ahh, now I can sleep.  Your analysis seems solid (although binomials were never my thing).
 
Jim

-----Original Message-----
From: insight-users-bounces at itk.org [mailto:insight-users-bounces at itk.org]On Behalf Of Stefan Klein
Sent: Thursday, March 17, 2005 8:27 AM
To: Insight-users at itk.org
Subject: RE: [Insight-users] random number generator


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  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/f54b6109/attachment.html