More QuickCheck and some geometric probability

by Anders Leino

Last time I didn’t write any customized generators, because the types I tested over were standard.
This time I will write two generators for unit vectors (aka norm vectors), i.e. vectors on the unit sphere.
There will also be some  pictures, so let’s include some OpenGL stuff as well.

> import Graphics.Rendering.OpenGL
> import Graphics.UI.GLUT
> import Test.QuickCheck
> import Test.QuickCheck.Gen
> import System.Random

I will keep it in three dimensions, for purposes of easily defined data-types and clarity.

> data Vector = Vector (Double, Double, Double)
>   deriving (Eq, Show)

A plain vector can be anywhere, so writing a generator is straight-forward.

> instance Arbitrary Vector where
>   arbitrary = do
>     arbitrary >>= return . Vector

Note that Vector takes a (Double, Double, Double).
QuickCheck provides instances of Arbitrary for Double and for Arbitrary a => (a,a,a), so there is not mutch to do but to package this up into a Vector.

Generating unit vectors

We are of course interested in the norm of a vector:

> norm (Vector (x,y,z)) = sqrt $ x**2 + y**2 + z**2

Now we can write more interesting generators, intended to generate vectors on the unit sphere.

Here is a norm generator which I call “cubical” for reasons which will soon be obvious.

> normGenCube = do
>   v@(Vector (x,y,z)) <- arbitrary `suchThat` (/= Vector (0.0, 0.0, 0.0))
>   return $ Vector (x/norm v, y/norm v, z/norm v)

This generator simply generates a non-zero vector as we defined above, and then scales it to have unit length.
Thus, one would expect the distribution to have cubical symmetry.

I initially expected this distribution to be denser near the projection of the corners of a cube onto the unit sphere, and less dense near the center of the projections of the faces of a cube.

However, the following picture tells the opposite story.
(The picture have been augmented by green points indicating the position of the unit cube face points on the sphere.)

> picture_1 = renderDistribution normGenCube (10^5) renderCubeFacepoints

(See below code if you want to find out how the rendering works.)

The distribution tends to bunch up around the face points instead, i.e. the very opposite of what I was expecting!

At the moment, I don’t exactly understand why this is.
To try and get to the bottom of it, I first went and looked at the QuickCheck manual.
Before long, I found myself digging trough the QuickCheck source code.

The QuickCheck source reveals that Arbitrary instances for triplets of Doubles are lifted in the most straight-forward way: we simply generate three numbers with arbitrarySizedFractional and put them into a triple. Nothing unexpected here.

The unexpected distribution must come from the way in which arbitrarySizedFractional is implemented.
To test this assumption, here is a normGenCube’ which uses the generator for Int instead of Double, and then simply converts it to a double using fromIntegral:

> normGenCube' = do
>   (x,y,z) <- arbitrary `suchThat` (/= (0, 0, 0)) :: Gen (Int, Int, Int)
>   let v@(Vector (x',y',z')) = Vector (fromIntegral x, fromIntegral y, fromIntegral z)
>   return $ Vector (x'/norm v, y'/norm v, z'/norm v)

And here is a picture.

> picture_2 = renderDistribution normGenCube' (10^5) renderCubeCorners

Nice! Avoiding the adjective ‘sized’ gives us the expected distribution.
I’ll leave the question of ‘Why?’ as an exercise to the reader.

Here is a different generator:

> normGenSphere = do
>   (theta,phi) <- arbitrary
>   return $ Vector (cos theta * cos phi, cos theta * sin phi, sin theta)

This generator is biased toward the north and south poles ((0,0,1) and (0,0,-1)) of the sphere.
(It is easy to understand why, by looking at the parametrization, so this generator behaves as we would predict.)

> picture_3 = renderDistribution normGenSphere (10^5) renderNorthSouthPoles

Making the pictures

Here is some quick and dirty OpenGL code which renders the pictures shown above.

For an nice introduction to OpenGL under haskell, see here, and also the sequel.

> renderDistribution gen num augment = do
>   getArgsAndInitialize
>   initialDisplayMode $= [RGBAMode, WithAlphaComponent]
>   createWindow "Plots"
>   blendFunc $= (SrcAlpha, OneMinusSrcAlpha)
>   blend $= Enabled
>   matrixMode $= Modelview 0
>   let eye = Vertex3 2.0 2.0 2.0
>       look = Vertex3 0.0 0.0 0.0
>       up = Vector3 0.0 0.0 1.0
>   lookAt eye look up
>   displayCallback $= do
>     clear [ColorBuffer]
>     renderSamples $ unGen (vectorOf num gen) (mkStdGen 12345) 1000
>     augment
>     flush
>   reshapeCallback $= (Just $ \s@(Size w h) -> do
>     viewport $= (Position 0 0, s)
>     matrixMode $= Projection
>     loadIdentity
>     perspective 90.0 (fromIntegral w / fromIntegral h) 0.0 10.0)
>   mainLoop
> -- render samples as semi-transparent grey points
> renderSamples vs = do
>   pointSize $= 1.0
>   color $ Color4 0.5 0.5 0.5 (0.2::Double)
>   renderPrimitive Points $ mapM_ (vertex . \(Vector (x,y,z)) -> Vertex3 x y z) vs
>   -- render helper points at north/south pole
> renderNorthSouthPoles = do
>   pointSize $= 3.0
>   color $ Color3 1.0 0.0 (0.0::Double)
>   renderPrimitive Points $ mapM_ vertex
>     [Vertex3 0.0 0.0 1.0, Vertex3 0.0 0.0 (-1.0::Double)]
> -- render points at projections of corners of
> -- unit cube onto unit sphere
> renderCubeCorners = do
>   pointSize $= 3.0
>   color $ Color3 1.0 0.0 (1.0::Double)
>   renderPrimitive Points $ mapM_ vertex
>     [Vertex3 (x/sqrt 3) (y/sqrt 3)  (z/sqrt 3)
>     |x <- [1.0::Double,-1.0], y <- [1.0,-1.0], z <- [1.0, -1.0]]
> -- render face-points at faces of unit cube
> renderCubeFacepoints = do
>   pointSize $= 3.0
>   color $ Color3 0.0 1.0 (0.0::Double)
>   renderPrimitive Points $ mapM_ vertex $
>     concat [[Vertex3 w 0.0 0.0, Vertex3 0.0 w 0.0, Vertex3 0.0 0.0 w] |
>             w <- [1.0,-1.0]::[Double]]
About these ads