

In [18]:

    
(require gamble
         racket/list)

The littlest radar blip problem.

This should illustrate the problems of representing uncertainty about existence, number, and origin.

I think a reasonable way to generate the blips at each time step is

generate the real blips, with each airplane having an independent, high probability of generating a blip, which has as its location a noisified version of the generating airplane's true location;
generate the number of false blips;
pick a random location for each false blip;
concatenate the true and false blips; and
randomly permute the order of the full blip list, to yield the observed blips.

We could achive this by randomly permuting the integers 0 to n-1, where n is the total number of blips, and using list-ref to get the right coordinate tuple.

To avoid possible type confusion, maybe the airplanes and blips should be data types, or at least gensyms.

A trickier variant could have new planes entering the box at random times.



In [ ]:



In [19]:

    
;; Worlds--sets of things, represented as integer indices.
(defmodel things
   (deflazy n-airplanes (poisson 5))
   (defmem (n-false-alarms time) (poisson 2))
   (deflazy airplanes (range n-airplanes))
   (defmem (false-alarms time) (range (n-false-alarms time))))



In [49]:

    
(open-model things)
(displayln airplanes)
(displayln (false-alarms 1))









    



(0 1 2 3 4 5 6)
(0 1)

Airplane dynamics

Initial locations (at time 0) are chosen uniformly in a (-1,1) box. Initial velocities are random normal.

Location at time $t$ is location at $t-1$ plus velocity at $t-1$ plus zero-mean noise.

Velocity at $t$ is velocity at $t-1$ plus zero-mean noise.



In [20]:

    
(defmodel dynamics
  (define location-noise-std-dev 0.05)
  (define velocity-noise-std-dev 0.1)

  (defmem (location plane dimension time)
          (if (= time 0) 
              (uniform -1 1)
              (normal (+ (location plane dimension (- time 1))
                         (velocity plane dimension (- time 1)))
                      location-noise-std-dev)))

  (defmem (velocity plane dimension time)
          (+ (normal 0 velocity-noise-std-dev)
              (if (= time -1) 
                  0
                  (velocity plane dimension (- time 1))))))

Observation model

This is tricky. We observe not only the locations of a number of blips, but also the fact there are exactly $n$ of them. Not sure how to get this in, yet.

I'm recalling some insight we had with "seismic", that once we sample the number of real observations, the number of false alarms is fixed by our data. What did we do there?



In [21]:

    
(defmodel observation
    
   (open-model things)
   (open-model dynamics)
          
   ;; On each dimension, the observed position will be the true position plus 
   ;; zero-mean noise.
   (define obs-noise-std-dev 0.01)

   ;; Each airplane may or may not generate a blip
   (define observation-prob 0.9)

   ;; Randomly select a subset of planes for observation.
   ;; NEEDS "airplanes" from "things" d-model.
   (defmem (observed-planes time)
           (filter (lambda (airplane) (flip observation-prob)) airplanes))
   
   ;; Randomly perturb the true locations of the observed planes to give 
   ;; observed locations.
   ;; NEEDS "location ..." from "dynamics" d-model.
   ;; BAD: this forces evaluation of *all* observations for the dim and time
   ;; (over all planes).  
   (defmem (observation-blip airplane dimension time)
      (normal (location airplane dimension time) obs-noise-std-dev))
   
   ;; Generate locations for false blips. 
   ;; NEEDS "false-alarms" from "things" d-model.
   (defmem (hallucination-blip h-number dimension time)
           (uniform -1 1))
   
   ;; Randomly permute indices to assign observations to real or hallucinated
   ;; sources
   ;; NEEDS "n-false-alarms" from "things" d-model.
   (defmem (blip-order time) 
           (sample (permutation-dist (+ (length (observed-planes time))
                                (n-false-alarms time)))))
          
   ;; A "soft" observation of the number of blips, so we don't throw out too many
   ;; rather good blip sets for not matching on exact number.
   ;; Ideally, the standard deviation would anneal to near zero.
   (defmem (about-how-many-blips time)
           (normal (vector-length (blip-order time)) 0.5))
          
            
   ;; Use blip order, observations, and hallucinations to define an observed-blip
   ;; function.  
   ;; Forces blip-order, which forces observed-planes and n-false-alarms, which can
   ;; lead to a blip number that is out of bounds.
   ;; ??!!! the act of observing the coordinates of a blip number is an implicit 
          ;;; observation that there are
          ;; at least that many blips!!!!
   (define (observed-blip blip-number dimension time)
     (let ([raw-idx (vector-ref (blip-order time) blip-number)])
       (if (< raw-idx (length (observed-planes time)))
           (observation-blip raw-idx dimension time)
           (let ([idx (- raw-idx (length (observed-planes time)))])
             (hallucination-blip idx dimension time)))))
          
    (define (blip-source-plane blip-number time)
       (let ([raw-idx (vector-ref (blip-order time) blip-number)])
           (if (< raw-idx (length (observed-planes time)))
               raw-idx
               "false alarm"))))

Let's make sure everything behaves as expected:

Ideally, we could test our beliefs not only about the what we expect some variables to look like, but also our beliefs about what variables have and do not have values, having asked for others. I think this means being able to peek in to the memo tables. I can imagine a standard display for rvs and indexed sets of rvs. Though textual output might be fine.



In [23]:

    
(open-model observation)

(printf "observed planes for step 3: ~a\n" (observed-planes 3))

(printf "\nobserved coordinates of observed planes at step 3:\n")
(map (lambda (pair) (printf "~a\n" pair))
        (map (lambda (n) 
               (list (observation-blip n 0 3 ) 
                     (observation-blip n 1 3))) (observed-planes 3)))

(printf "\nobserved coordinates of hallucinated planes at step 3:\n")
(map (lambda (pair) (printf "~a\n" pair))
        (map (lambda (n) 
               (list (hallucination-blip n 0 3 ) 
                     (hallucination-blip n 1 3))) (false-alarms 3)))

(printf "\nblip-order for step 3: ~a\n" (blip-order 3))

(printf "\nabout how many blips for step 3: ~a\n" (about-how-many-blips 3))

(printf "\nobserved blip locations at step 3:\n")
(map (lambda (pair) (printf "~a\n" pair))
        (map (lambda (n) (list (observed-blip n 0 3 ) (observed-blip n 1 3)))
             (range (vector-length (blip-order 3)))))

(printf "\nblip sources for blips at step 3:\n")
(map (lambda (blip-n)
       (printf "~a\n" (blip-source-plane blip-n 3)))
     (range (vector-length (blip-order 3))))









    



observed planes for step 3: (0 2 3 4 5 6)

observed coordinates of observed planes at step 3:
(0.3927818047097665 0.7995508900442223)
(0.43072672082268 -0.9303208890547434)
(-0.37592817990791705 -0.4322011979362125)
(-0.7775057356388749 0.6381384083952282)
(-0.49647579038875905 -0.7002570495528102)
(-0.8354163917597375 0.5179818294481732)

observed coordinates of hallucinated planes at step 3:
(-0.30521747411350597 0.5594072026099773)
(0.6541575212182394 0.9157708735392294)

blip-order for step 3: #(3 5 7 4 0 1 6 2)

about how many blips for step 3: 8.30035226792435

observed blip locations at step 3:
(-0.37592817990791705 -0.4322011979362125)
(-0.49647579038875905 -0.7002570495528102)
(0.6541575212182394 0.9157708735392294)
(-0.7775057356388749 0.6381384083952282)
(0.3927818047097665 0.7995508900442223)
(-0.3637001156904466 0.3523995445515496)
(-0.30521747411350597 0.5594072026099773)
(0.43072672082268 -0.9303208890547434)

blip sources for blips at step 3:
3
5
false alarm
4
0
1
false alarm
2






    Out[23]:





(#<void> #<void> #<void> #<void> #<void> #<void> #<void> #<void>)

Some data

Let's have blip data in a table (list of lists) with columns "time", "blip number", "x", and "y".

Four real planes are in play:

Starts at the top left ((-1,1)) and moves right and slightly down.
Starts at the origin and moves right and up.
Starts at (1,0) and moves left and down.
Starts at (1,-1) and moves left and up.

Zero to three hallucinations are thrown in on each time step.

A function in the sampler makes an observation about the number of blips at each time, and a pair of observations for each blip's x and y coordinates at each time.



In [22]:

    
(define blip-data
  '((1 2 -1.0 0.99)
    (1 1 -0.01 -0.02)
    (1 3 0.99 0.001)
    (1 5 0.985 -0.99)
    (1 4 -0.9 0.2)
    
    (2 5 -0.96 0.97)
    (2 2 0.01 0.03)
    (2 6 0.96 -0.01)
    (2 3 0.97 -0.95)
    (2 4 -0.4 0.86)
    (2 1 0.2 0.67)
    
    (3 4 -0.9 0.91)
    (3 2 0.05 0.08)
    (3 3 0.91 -0.055)
    (3 1 0.92 -0.9)
    (3 5 0.2 -0.6)
    
    (4 6 -0.84 0.87)
    (4 3 0.09 0.11)
    (4 1 0.88 -0.15)
    (4 5 0.85 -0.84)
    (4 4 -0.2 0.99)
    (4 2 -0.5 -0.9)
    
    (5 1 -0.8 0.834)
    (5 2 0.13 0.17)
    (5 3 0.84 -0.19)
    (5 4 0.8 -0.8)
    (5 5 0.3 0.7)
    
    ))

(define times (sort (remove-duplicates (map first blip-data)) <))


(define (blips-for t)
  (filter (lambda (b) (= (first b) t)) blip-data))

(define (n-blips t)
  (length (blips-for t)))

(define blip-n (lambda (blp) (- (second blp) 1)))
(define x-val third)
(define y-val fourth)

The sampler

What would we like to find out, given our blip sequence?

Which blips are hallucinations (probably) and which are real?
Which blips from different times are generated by the same object? [hard]

Because of symmetries (...more), it won't do to query for which plane blips are generated by. It is more robust (if more costly and verbose) to ask which blips in consecutive time steps are generated by the same object (whichever that is). We're going to need a function to get the source object for a blip. [done--check]

We could show a sequence of plots, each of which shows

The blips recieved at that time, each color-coded for probability of true (blue) vs. hallucinated (red). For this, we need to save out the source of each blip, on each sample.
The estimated position of each numbered aircraft--just a cloud of dots, for now, since the numbers are not identifiable. For this, we need to save out a list (of variable length) of all the airplane locations, on each sample. We aren't averaging these, so the variable length should not be a problem.



In [23]:

    
(define radar-sampler
  (mh-sampler
   
   (open-model observation) ; this opens "things" and "dynamics"
   
   (for ([t times])
     (begin
       (observe (+ (length (observed-planes t)) (n-false-alarms t)) (n-blips t))
       (for ([blp (blips-for t)])
         (observe (observed-blip (blip-n blp) 0 t) (x-val blp))
         (observe (observed-blip (blip-n blp) 1 t) (y-val blp)))))

   
   (location 0 0 1)))



In [24]:

    
(generate-samples radar-sampler 10 #:burn 10 #:thin 10)









    Out[24]:





#(-0.17465999243333902 -0.17465999243333902 -0.17465999243333902 -0.17465999243333902 -0.17465999243333902 -0.17465999243333902 -0.17465999243333902 -0.17465999243333902 -0.17465999243333902 -0.16636515198044283)



In [ ]:

    
;; junk
   (defmem (source-plane blip)
           (let ([p-real (/ n-airplanes (+ n-airplanes (n-false-alarms))])
             (if (flip p-real) (discrete n-airplanes) "false alarm")))
   (defmem (blip-location blip dimension)
           (let ([s (source-plane blip)])
             (if (string=? s "false alarm")
                 (uniform -1 1)
                 (normal (location s dimension (time-of blip)) 0.01))))
                         
                         
                         
   (observe (+ (length (observed-planes 1)) (n-false-alarms 1)) 5)
   
   (observe (observed-blip 0 0 1) 0.5)
   (observe (observed-blip 0 1 1) 0.7)
   (observe (observed-blip 1 0 1) -0.5)
   (observe (observed-blip 1 1 1) -0.2)
   (observe (observed-blip 2 0 1) 0.15)
   (observe (observed-blip 2 1 1) 0.79)
   (observe (observed-blip 3 0 1) -0.95)
   (observe (observed-blip 3 1 1) -0.02)
   
   (observe (+ (length (observed-planes 2)) (n-false-alarms 2)) 5)
   ;(observe (about-how-many-blips 2) 7)
   (observe (observed-blip 0 0 2) 0.53)
   (observe (observed-blip 0 1 2) 0.67)
   (observe (observed-blip 3 0 2) -0.45)
   (observe (observed-blip 3 1 2) -0.22)
   (observe (observed-blip 1 0 2) 0.18)
   (observe (observed-blip 1 1 2) 0.89)
   (observe (observed-blip 2 0 2) -0.85)
   (observe (observed-blip 2 1 2) -0.11)



In [31]:

    
(sort (remove-duplicates '(5 2 3 3 5 2 3)) <)









    Out[31]:





(2 3 5)



In [ ]: