The Coq mailing list

[AUGER Cedric <Cedric.Auger AT lri.fr>]

Le Tue, 24 May 2011 10:50:44 +0200,
Fontaine Allyx 
<fontaine AT labri.fr>
 a écrit :

> Le 23/05/2011 11:18, AUGER Cedric a écrit :
> > Le Sun, 22 May 2011 17:41:50 +0200,
> > David 
> > Baelde<david.baelde AT gmail.com>
> >   a écrit :
> >
> >> Hi Allyx,
> >>
> >> If you haven't already, you should have a look at the paper by
> >> Audebaud&  Paulin (http://dx.doi.org/10.1016/j.scico.2007.09.002)
> >> about the ALEA library. It starts with a gentle introduction to the
> >> monad that is used, detailing its advantages and drawbacks.
> >>
...
> >>
> >> Long story short, Pierre Corbineau's solution will work. However, I
> >> wouldn't bother to make extraction work for such small examples:
> >> just rewrite the code in ML yourself.
> >>
> >> Have fun,
> >
> > I also guess that if "'a" is a type which cardinality is quite
> > small, and you can extract [0,1] to float, then you can write a
> > short ML code:
> >
> > let dice_throw =
> >    let sieve = Random.float 1. in
> >    let rec get_rand acc l =
> >      match l with
> >      | v::l ->  let acc = acc + μ v in
> >                if acc>  sieve then v else get_rand acc l
> >      | [] ->  failwith "Hmm, this is unexpected!"
> >    in get_rand 0 [list of all the elements of your type]
> >
> > of course if you have big types for 'a such as int it can be too
> > long and even false due to rounding errors (the failwith case can
> > even occur).
> > You also have to identify the "μ" function (if it is the one I think
> > of, it is also called μ in Alea, so it should be extracted to its
> > Unicode code, as unicode is not supported by OCaml, leading to some
> > awfull identifier).
> > Even in this case, I am not sure it can be very usefull, but if it
> > is, then you won't have to rewrite all the ML code, and won't have
> > to play with the Extraction mechanism too much (and using P.
> > Corbineau solution seems also hazardous, as you change the
> > extracted objects to objects which are semantically very different;
> > furthermore, using side effects can mess a lot with extraction if
> > some optimizations are performed ("let x = side_effect () in let y
> > = side_effect () in ..." could be extracted to "let x = side_effect
> > () in let y = x in ..." as well as identity due to Coq semantics,
> > but the two functions are clearly different with OCaml)).
> 
> Thanks for your answers!
> 
> Actually, as I have several examples more or less simple, I would
> like to have an automatic way which ensures that Ocaml functions
> really corresponds to Coq functions. But as I don't know well all
> those mechanisms, I indeed begin with writting ML code by myself and
> having a look at extraction.
> 
> I am not sure of my understanding of the role of the accumulator in 
> dice_throw. Is it needed to return an element of my type uniformly at 
> random thanks to the μ function? Does the failwith case not be
> reached because the sum of the probabilities is supposed to be
> greater than 1?
> 
> Allyx

Almost yes for the last paragraph (in fact the sum is expected to be
exactly 1, and not "greater" than 1; in french "greater" has a large
meaning, but in english it has a strict one).

In fact I never used the Alea library, so what I wrote may be wrong,
but if you want some commented version of the function I wrote:
(with the notation:
Int(a,b,P) for:
 b
∫ P(x) ∂x
 a
and:
f^(-1) for:
f⁻¹
)

 let dice_throw =
    (*First we throw a dice on [0,1] with uniform probability,
      assuming [Random.float] is able to do so,
      which is not the case in fact as it is a pseudo-random
      generator and float doesn't cover all reals in [0,1],
      but it is expected to be a good approximation of this law.
      Let us bind this result to the "sieve" identifier.
    *)
    let sieve = Random.float 1. in
    (* Now thanks to theorems which I do not remember the names,
       we have some transfer function between uniform law and
       any Riemann integrable law P on [a, b] such that Int(a,b,P)=1
       (I don't recall if it can be generalized to Lebesgue,
        but that doesn't matter here) with:
        X |--> (x |--> Int(a, x, P))^(-1)(X)
       So that now we want to compute "(x |--> Int(a, x, P))^(-1)"
    *)
    let rec get_rand acc l =
      (* [acc] is meant to be "Int(a, predecessor(List.hd l), P)"
         where [a] is the first element of the type
         (in fact as the type is finite, the integral is just
          a finite sum)
      *)
      match l with
      | v::l ->  let acc = acc + μ v in (*μ is our P law*)
                if acc>  sieve
                then v (*ok, we have reversed our integral*)
                else get_rand acc l
                 (*we didn't reverse our integral and had to go
                   further*)
      | [] -> (*if this point is reached, then by definition
                of probability laws, we should have that [acc]=1.
                and so, due to definition of sieve we should have
                halted earlier. Of course due to rounding,
                we could have [acc]=0.99988665 or so and
                [sieve] = 0.9999039 or so, so that this case
                may happen if we are very unlucky
              *)
              failwith "Hmm, this is unexpected!"
   in
   get_rand 0 [list of all the elements of your type]


Hope these comments are not too messy!
Once again, I guess it is what you would like, but I am not sure it is
very convenient. And if you want to do big benchmarks to be able to
make some observations, it may be quite slow. But I have no really
better solution.