sampling a pmf

Stream: maths

[Ben Ryjikov]

Hi everyone, is there a way to sample a value from a probability mass function without creating a new pmf?

I'm using mathlib's pmf definition. The one way I've found to sample is pmf.bind, which samples from a pmf, but also passes the sampled value into a function that creates a new pmf. I would like to sample a value without creating a new pmf, which means bind isn't exactly what I'm looking for.

One solution I've found is to have bind sample a value, then construct the pmf which always returns that value (constructed by pmf.pure). However, I then would like to extract the value from that pmf, which I haven't found a way to do.

Can this be done? If so, how?

(also if this has been discussed or should be in a different stream, please redirect me there (I'm a new member please bear with me))

[Adam Topaz]

the docs#pmf.bind is the monadic bind. As is usual with monads, there is no way to "extract" the values, but you can "extract" them if you stay within the monad environment.

[Adam Topaz]

If you have a random number generator and some "computable" pmf on X, then in principle it is possible to "sample" from such a pmf, but the result would be something of type io X, not X. But I don't know if anything like this has been implemented.

[Adam Topaz]

import probability.probability_mass_function.monad

noncomputable
example (P Q : pmf ℕ) : pmf ℕ := do {
  p ← P, -- "sample" from P
  q ← Q, -- "sample" from Q
  return (p + q) -- return the sum of the samples.
}

[Eric Wieser]

I think even if you build your own pmf from scratch you'd still have a really hard time making it computable because the monad framework that enables do notation doesn't let you put constraints on your types.

[Ben Ryjikov]

Ah I see. Thanks for the feedback.

[Adam Topaz]

src#pmf

[Alex J. Best]

Eric Wieser said:

I think even if you build your own pmf from scratch you'd still have a really hard time making it (edit: the monad instance) computable because the monad framework that enables do notation doesn't let you put constraints on your types.

What constraints are you thinking of? I guess I don't quite get what you're saying here

[Adam Topaz]

Here's some fun with Lean4:

inductive PreGiry : Type _ → Type _ where
  | uniform : Nat → Nat → PreGiry Nat
  | bind : PreGiry α → (α → PreGiry β) → PreGiry β
  | dirac : α → PreGiry α

instance : Monad PreGiry where
  pure := PreGiry.dirac
  bind := PreGiry.bind

def PreGiry.sample : PreGiry α → IO α
  | uniform m n => IO.rand m n
  | bind P f => P.sample >>= fun t => (f t).sample
  | dirac x => return x

-- TODO: Change 2000000 to some other value?
def Uniform : PreGiry Float := do
  let p ← PreGiry.uniform 0 2000000
  return p.toFloat / 2000000

-- Probably good enough ;)
def Pi := 3.141592653589793238462643383279502884

-- The Box-Muller transform
def Normal : PreGiry Float := do
  let u1 ← Uniform
  let u2 ← Uniform
  return Float.sqrt (-2 * Float.log u1) * Float.cos (2 * Pi * u2)

#eval Normal.sample

[Scott Morrison]

Fun. Could you explain why you called it PreGiry?

[Eric Wieser]

What constraints are you thinking of?

The example that comes to mind is making docs#finset.bUnion a has_bind instance. It's not possible because you can't say "finset is a monad over types with decidable equality". This thread has more details.

[Adam Topaz]

I named it PreGiry because there is a map of (possibly non-lawful) monads from this to the Giry monad. This map has a dense image (in some sense that I don't know how to make precise right now) and is (very) noninjective. Thus in some sense the Giry monad is some completion of some quotient of PreGiry. Making all this precise could be a fun (probably nontrivial) exercise!