F# で最高事後密度区間 (highest posterior density interval, HPDI) を求めます。

let hpdi mass (values : float list) =
   let size = List.length values
   let shift = int <| mass * float size
   let sorted = List.sort values
   let width = Seq.map2 (-) (Seq.skip shift sorted) sorted |> Seq.toList
   let minIndex =
      Seq.zip (Seq.initInfinite id) width
      |> Seq.minBy snd
      |> fst
   sorted.[minIndex], sorted.[minIndex + shift]

取りうる区間の中でもっとも幅が狭いものを選択しています。

例として標準正規乱数を 100,000 個生成してその 95% HPDI を求めてみます。

let rec normal (random : System.Random) =
   seq {
      let u1 = random.NextDouble ()
      let u2 = random.NextDouble ()
      let r = sqrt <| -2.0 * log u1
      let theta = 2.0 * System.Math.PI * u2
      yield r * sin theta
      yield r * cos theta
      yield! normal random
   }
let samples = normal (System.Random ()) |> Seq.take 100000 |> Seq.toList
hpdi 0.95 samples |> printfn "%A"  // (-1.940508623, 1.96964806)

結果がおよそ (-1.96, 1.96) になっていますね。