Draw-bability!

I started playing with Typed Racket, which is basically an attempt to build a type system that can accommodate idiomatic Scheme programming without modification. These leads them to such novel features as untagged unions, non-uniformly typed varargs, and (pretty helpful when you have untagged unions),  a type-refining version of instanceof called “occurrence typing”.

As more of an F# programmer than a Scheme programmer, the type system feature that most clearly appealed to me was the ability to abstract over variable numbers of function arguments. You see functions of the form foo2, foo3, foo4,… all the time, identical except for their number of arguments – a clear failure of the type system’s ability to abstract out boilerplate. Typed Racket then gives us an exciting opportunity to explore relatively uncharted territory, modes of abstraction that are not possible in mainstream functional programming languages. One possible use of this feature is to create some really clean foundations for a probabilistic modeling language. I’ve found that working with the data type of probability distributions really showcases the power and expressivity of functional programming, and it’s going to be fun to see if this type system can help us come up with an even cleaner formulation.

Here’s our definition of a probability distribution and a couple helper functions:

(define-type (Dist a) (-> a))
(: sample-dist (All (a) ((Dist a) -> a)))
(define (sample-dist d) (d))
(: always (All (a) (a -> (Dist a))))
(define (always val) (λ () val))

Here, probability distributions of type a are represented as non-referentially-transparent “sampling functions” that take no arguments, and return a sample value of type a. The usual way to make these into pure functions is to change the type from () -> a to  [0, 1] -> a, and feed them a random source, but the version above is even simpler and fine for our purposes.

(always some-value)

produces a probability distribution which is always some-value, that is, it returns some-value every time you sample it.

The type Dist is a functor (among many other things), so it has an operation, fmap:

(: dist-fmap (All (a b) ((a -> b) -> (Dist a) -> (Dist b))))
(define ((dist-fmap func) d)
  (λ () (func (sample-dist d))))

fmap is great for probability modeling – it lets you “lift” operations up into the world of distributions. You can take a function of normal values and turn it into a function over random values. You can see the analogy to lists and mapmap takes a normal function and turns it into a function from lists to lists. But our dist-fmap only works for functions of one argument. Boo. So now we can’t even do things like add two random numbers and get the resulting distribution. We could go ahead and define a bunch more functions, fmap2, fmap3, etc, that let us lift functions of varying numbers of arguments into the world of distributions. But Typed Racket has a trick up its sleeve:

(: dist-fmap* (All (c a b ...) ((a b ... b -> c) ->
                                (Dist a) (Dist b) ... b -> (Dist c))))
(define (dist-fmap* f)
  (λ (da . dbs) (λ () (apply f (sample-dist da) (map sample-dist dbs)))))

This definition uses a “dotted type” to represent the type of all fmap-functions for distributions. The type variable b ... b represents a list of type variables, each of which can be substituted with any type, forming a truly general definition of fmap. Of course this means Dist isn’t just a functor any more, it has a bit more structure, which is a subject for a blog post in and of itself. I like to give this function a shorthand:

(define ^ dist-fmap*)

Which to me conveys the “lifting” idea of the operator. With the power of the magic ^, we can now take any function and turn it into a function on distributions. But wait, we don’t want everything to be a distribution, sometimes we need some constant values. While we think of ^ as allowing us to lift functions into the land of probability, we can think of always as the ability to lift constant values. So one more shorthand:

(define p always)

and now we can easily lift functions and values (I know, functions are values, but you know). Well, now that we have these cool tools, let’s see what they can do. Let’s start with the uniform random variable on the interval [0,1]:

(: uniform (Dist Real))
(define uniform random)
> (uniform)
0.5597565477782307
> (uniform)
0.7969202147236572

Let’s kick it up a notch:

> (((^ min) uniform uniform uniform uniform uniform)))
0.06934773489468006
> (((^ min) uniform uniform uniform uniform uniform)))
0.09233879582175743
> (((^ min) uniform uniform uniform uniform uniform)))
0.14803391736723828

Cool! That’s called the 5th order-statistic. The nth order-statistic is the distribution of the smallest element in n samples of a uniform random variable. And it’s written in about the cleanest way you could imagine. But repetition is the mother of crappy code, and having to repeat “uniform” all those times is lame. That gives me an idea:

(: lift-list (All (a) ((Listof (Dist a)) -> (Dist (Listof a)))))
(define (lift-list ds)
  (foldl (^ cons) (p '()) ds))

(: sample-population (All (a) ((Dist a) Natural -> (Dist (Listof a)))))
(define (sample-population d sample-size)
  (lift-list (repeat d sample-size))) 

(: nth-order-statistic (Natural -> (Dist Real)))
(define (nth-order-statistic n)
  (: min-list ((Listof Real) -> Real))
  (define (min-list xs) (apply min xs))
  ((^ min-list) (sample-population uniform n))

We could have solved this problem a bunch of ways, but lift-list and sample-population are useful for more than just this one example. They allow us to create a derived distribution of sample sets that we can then fold in interesting ways. min is one of those folds, and the combination of sample-population and lifted min gives us the the order statistics. We can even calculate sample mean and standard deviation functions, parametrized by sample size (bodies of statistical functions mean and sigma omitted for brevity):

> (((^ mean) (sample-population d 1000))))
0.5139713415117088
> (((^ mean) (sample-population d 1000))))
0.512734382343886
> (((^ sigma) (sample-population d 1000))))
0.29030493445369654
> (((^ sigma) (sample-population d 1000))))
0.28459856532796873

Seems pretty accurate (and totally meta – distributions of statistics!).

Anyhow, we’ve established that we can do some pretty cool numerical modeling with this functional language, I think. So far though, this is pretty tame. With the ^ operator I’m making fun use of TR’s typed variadic arguments, but I want to get a bit more adventurous and a bit more Rackety. I really enjoyed the Pictures tutorial, so I wanted to do something in that vein. Unfortunately, slideshow, the library used, is not yet ported to Typed Racket. I would love to port it myself, but I was having trouble digging around and figuring out everything I had to wrap up, and dammit, I want to see some pictures, so let’s take a shortcut. Luckily we have the magic typed/racket/no-check switch, which is going to let us take this existing typeful infrastructure we have created, and link in some dynamically typed code. This should be a great existence proof of the fact that programming in a dynamic language does not mean programming without types. So, we put the following lines at the top of the page:

#lang typed/racket/no-check
(require slideshow)

And we’re off to the races. The types of all the following expressions we’re going to write are (Dist Pict), where Pict is a picture data type. That’s right, we’re dealing in picture-valued random distributions now. Once again for the record, although we’re not officially type-checking, we know exactly what types we’re dealing with and are writing perfectly type-safe code. Here’s a circle of radius 0-50, evenly distributed:

> (((^ circle) ((^ *) (p 50) uniform)))

> (((^ circle) ((^ *) (p 50) uniform)))

> (((^ circle) ((^ *) (p 50) uniform)))

> (((^ circle) ((^ *) (p 50) uniform)))

> (((^ circle) ((^ *) (p 50) uniform)))

Now maybe I’m just sheltered, but this totally rocks my world. Here we’ve created this beautiful little mini-language for creating probabilistic data, and we were in some ways guided there by the types. The design feels really clean and sort of inevitable (modulo minor quibbles like restoring referential transparency, etc.) – I mean to me this feels as about as good a way of drawing random pictures as I could have imagined.

Anyhow, let’s draw some more pictures. You get some pretty cool patterns by overlaying bunches of these circles:

> ((apply (^ cc-superimpose)
          (repeat ((^ circle) ((^ *) (p 100) uniform)) 10)))

> ((apply (^ cc-superimpose)
          (repeat ((^ circle) ((^ *) (p 100) uniform)) 10)))

> ((apply (^ cc-superimpose)
          (repeat ((^ circle) ((^ *) (p 100) uniform)) 10)))

Neato. I have to say, I think this gives me an opportunity to combine functionality from different libraries in ways the authors truly hadn’t anticipated. Let’s try those drawings again, but with the Gaussian distribution imported from the “random” library on PLaneT:

(require (planet schematics/random:1:0/random))
> ((apply (^ cc-superimpose)
          (repeat ((^ circle) ((^ *) (p 60) ((^ abs) random-gaussian))) 20)))

> ((apply (^ cc-superimpose)
          (repeat ((^ circle) ((^ *) (p 60) ((^ abs) random-gaussian))) 20)))

> ((apply (^ cc-superimpose)
          (repeat ((^ circle) ((^ *) (p 60) ((^ abs) random-gaussian))) 20)))

>

got to make sure we take the absolute value (the lifted absolute value, (^ abs), naturally), as negative circle radii are a big turn-off for drawing libraries. So, it’s easy to see that with the Gaussian, the radii are much more concentrated around the mean.

Anyhow, I’m going to stop for now with the vaguely suggestive drawings and take a breather. Look out for another post getting into the fmap*/^ operator some more and talking about why I think its cool. Also, more pictures. This stuff is pretty weak sauce compared to what could be done – the libraries are very powerful. I’m also quite interested in trying to draw some less geometric things. Imagine drawing a bunch of stick figures, eyes, mouths, and limbs at constrained random angles, and seeing what kind of scenes emerged? Stochastic cartooning, here we come! I’m sure I’m not the first to think of this sort of thing – does anyone know of similar projects?

Advertisements

One Response to Draw-bability!

  1. […] post on using HLists to kinda-sorta simulate variable-arity polymorphism, as a followup to my very frist psot using Typed Racket’s variable-arity features to make a “fmap*” function that […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: