haskell - 在 Haskell 中实现 Metropolis 算法

Question

我正在尝试编写一个基于 Haskell 中的 Metropolis 算法的 MCMC 程序，但我在从概率分布中采样(生成伪随机数)和构建程序时遇到问题。目前，我很高兴使用带有硬编码种子的生成器，而不是处理 IO 的复杂性。

看来我应该使用状态单子(monad)来跟踪随机生成器状态、先前的马尔可夫链状态、卡方值以及算法每一步之间的接受计数，然后最终收集所有马尔可夫链状态和最后接受计数。这是最好的/惯用的方法吗？如果是这样，程序的布局应该是什么(即提案函数和都会步骤函数的类型签名等)。

我见过一些处理随机数的示例程序，其中特定长度的随机数列表是从某种概率单子(monad)生成的，然后通过一些简单的函数来执行计算。如果可能的话，我真的很想避免这种由内而外的程序形式。

编辑:暂时删除 WIP 代码。

Answer 1

这里是一些关于编写惯用的 Haskell 的反馈。

• 除非您正在编写单子(monad)代码，否则请在纯函数中使用 do(即constructMuTable、metropolis)是非常不惯用的Haskell。

  而不是

  foo = do
      let x = ...
          y = ...
          z = ...
      bar x y z
  

  直接写

  foo =
      let x = ...
          y = ...
          z = ...
      in bar x y z
  

  或使用where而不是let ... in ...。

• eta 减少。在某些地方(main 中的 zVec、muVec、sigmaVec)，您编写了 (\x -> f x)。这相当于 f、模 _|_、seq 等。

• 使用Data.Vector.Unboxed。您有很多 V.Vector Double，它存储盒装 Double，效率可能很低。对于像 Double 这样的基本类型，使用未装箱的向量(可能)使用更少的内存来获得更快的代码。

• 如果可能，请避免使用 (!!) 索引列表。使用 Data.Vector 代替，因为 V.! 是 O(1)，而 (!!) 是 O(n)。

• 看来您可以在此处使用 State monad 来清理代码。然而，以目前的命令式形式，我很难看到这种转变。

也许您可以尝试应用我给出的一些建议并简化一些大型、密集的函数，然后对您的算法的更高级别的反馈将变得更加明显。