我如何在Julia中编写任意连续分布，或者至少模拟从一个连续分布中抽样?

Question

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

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

我如何从一个随机变量的样本值与此PDF在朱莉娅?(或者，我如何至少从这样一个随机变量中模拟抽样?)

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

(这已经是离散分布的答案了，我如何在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]区间内的分布。为了回答这个问题，我假设你的概率分布函数在每个区间上都是均匀的。鉴于此，我认为实现您自己的发行版的最简便的方法是创建一个新的子类型，在本例中是ContinuousUnivariateDistribution。示例代码如下:

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随机数的最小值。查看这里的其他方法列表，您可能想要添加到您的新发行类型(取决于您的用例，也许您不会麻烦)。

注意，如果效率是一个问题，你可以研究一些方法，这些方法将允许你最小化d.theta^2操作的数量，例如Distributions.sampler。或者，您可以将theta^2内部存储在MyDistribution中，但始终显示底层theta。真的由你决定。

最后，在函数输出上并不需要类型注释。我只是为了清楚地把它们包括进来。