julia - 我如何在 Julia 中编写任意连续分布，或者至少模拟一个连续分布的采样？

Question

假设我有一个任意概率分布函数 (PDF)，定义为函数 f，例如:

using Random, Distributions

#// PDF with parameter θ ϵ (0,1), uniform over 
#// segments [-1,0) and [0,1], zero elsewhere
f = θ -> x -> 
  (-1 <= x < 0) ? θ^2 :
  (0 <= x <= 1) ? 1-θ^2 :
  0

如何在 Julia 中使用此 PDF 从随机变量中采样值？ (或者，我至少如何模拟从这样的随机变量中进行采样？)

即我想要相当于 rand(Normal(),10) 的(标准)正态分布中的 10 个值，但我想使用函数 f 来定义分布使用(类似于 rand(f(0.4),10) - 但这不起作用)

(这已经是How can I write an arbitrary discrete distribution in Julia?处离散分布的答案:但是我想使用连续分布。在https://juliastats.org/Distributions.jl/v0.14/extends.html处创建采样器有一些细节，我认为这可能有用，但我不知道不明白如何应用它们。同样在 R 中，我使用了 https://blogs.sas.com/content/iml/2013/07/22/the-inverse-cdf-method.html 中描述的逆 CDF 技术来模拟此类随机变量，但我不确定如何在 Julia 中最好地实现它。)

Answer 1

第一个问题是，您提供的并不是概率分布的完整规范，因为它没有说明区间 [-1, 0) 内或区间内的分布。间隔[0, 1]。因此，出于本答案的目的，我将假设您的概率分布函数在每个间隔上都是均匀的。鉴于此，我认为实现您自己的分布的最朱利安方法是创建一个新的子类型，在本例中为 ContinouslyUnivariateDistribution。示例代码如下:

using Distributions
struct MyDistribution <: ContinuousUnivariateDistribution
    theta::Float64
    function MyDistribution(theta::Float64)
        !(0 <= theta <= 1) && error("Invalid theta: $(theta)")
        new(theta)
    end
end
function Distributions.rand(d::MyDistribution)::Float64
    if rand() < d.theta^2
        x = rand() - 1
    else
        x = rand()
    end
    return x
end
function Distributions.quantile(d::MyDistribution, p::Real)::Float64
    !(0 <= p <= 1) && error("Invalid probability input: $(p)")
    if p < d.theta^2
        x = -1.0 + (p / d.theta^2)
    else
        x = (p - d.theta^2) / (1 - d.theta^2)
    end
    return x
end

在上面的代码中，我为新分布实现了 rand 和 quantile 方法，这是能够进行像 rand( MyDistribution(0.4), 20) 从新分布中抽取 20 个随机数。请参阅here有关您可能想要添加到新的分发类型中的其他方法的列表(取决于您的用例，也许您不会打扰)。

请注意，如果效率是一个问题，您可以研究一些可以最大程度地减少 d.theta^2 操作数量的方法，例如Distributions.sampler。或者，您可以将 theta^2 存储在 MyDistribution 内部，但始终显示底层 theta。真的取决于你。

最后，您实际上并不需要函数输出上的类型注释。为了清楚起见，我只是将它们包括在内。