用Haskell实现Metropolis算法

Question

我正在尝试用Haskell编写一个基于Metropolis算法的MCMC程序，我遇到了从概率分布中采样(生成伪随机数)和组织程序的问题。就目前而言，我很乐意使用具有硬编码种子的生成器，而不是处理复杂的IO。

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

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

编辑:暂时删除在制品代码。

Answer 1

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

• ，除非您正在编写一元代码，在纯函数中使用do (即。constructMuTable，metropolis)非常不地道。

而不是

foo = do，让x= ...Y= ...Z= ...条形x y z

只需写

foo =让x= ...Y= ...Z= ...在条x y z中

或者用let ... in ....

• eta-reduce.代替where在某些地方(main中的zVec、muVec、sigmaVec )，您已经编写了(\x -> f x)。这等同于f，modulo _|_，seq等。

• 使用Data.Vector.Unboxed。您有很多存储盒装Doubles的V.Vector Double，效率可能很低。对于像Double这样的原语类型，如果可能的话，在(!!)中使用更少的memory.

• avoid索引表，为(潜在的)更快的代码使用未装箱的向量。使用Data.Vector代替，因为V.!是O(1)，而(!!)是O(1)，似乎您可以在这里使用State monad来清理您的代码。然而，在目前的命令式形式下，我很难看到这种转变。

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