从 Julia 中的自定义概率密度函数中抽取随机数

Question

我想从修改后的指数分布中抽取随机数：

p(x) = C * a * Exp[-(a*x)^b] 并使用 C=1/Gamma[1 + 1/b] 进行归一化。

我怎样才能在朱莉娅中做到这一点？不幸的是，我对 Julia 的经验很少，也没有创建自定义随机数的经验。我将非常感谢任何帮助。

Answer 1

如果我没记错的话，这是一个p-广义高斯分布，它在Distributions.jl中有一个相当有效的实现：

using Distributions
mu = 0       # your location parameter
alpha = 1/a  # your scale parameter
beta = b     # your shape parameter
p = PGeneralizedGaussian(mu, alpha, beta)

使用单变量分布的 Distributions.jl API，您可以通过将其传递给导出rand方法来从此分布中进行采样。PGeneralizedGaussian下面是如何使用mu = 0、alpha = 1/2和 来从分布中采样五个独立标量的示例beta = 3：

julia> p = PGeneralizedGaussian(0, 1/2, 3);

julia> rand(p, 5)
5-element Vector{Float64}:
  0.2835117212764108
 -0.023160728370422268
  0.3075395764050027
 -0.19233721955795835
  0.21256694763885342

如果您想尝试自己实现分布（我不建议这样做，除非您将其作为 Julia 中数学编程的练习），您需要定义一个类型来保存分布的静态参数（在您的情况下，只需形状参数和比例参数），然后定义并导出此处列出的方法，以将 Distributions.jl API 扩展到您的自定义分布。特别是，您需要定义：

struct MyDistribution <: ContinuousUnivariateDistribution
  # ... your distribution parameters here
end

rand(::AbstractRNG, d::MyDistribution) # sample a value from d
sampler(d::MyDistribution)             # return d or something that can sample from d more efficiently
logpdf(d::MyDistribution, x::Real)     # compute the log of the pdf of d at x
cdf(d::MyDistribution, x::Real)        # compute the cdf of d at x
quantile(d::MyDistribution, q::Real)   # compute the qth quantile of d
minimum(d::MyDistribution)             # return the minimum value supported by d
maximum(d::MyDistribution)             # return the maximum value supported by d
insupport(d::MyDistribution, x::Real)  # query whether x is supported by d

该包的文档非常好，因此如果您想学习 Julia，这是一个很好的入门方法。