从Julia中的自定义概率密度函数中生成随机数。

huangapple go评论59阅读模式
英文:

Draw random numbers from a custom probability density function in Julia

问题

我想从修改后的指数分布中随机生成数字:

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

我该如何在Julia中实现这个?不幸的是,我对Julia的经验有限,没有创建自定义随机数的经验。我将非常感激任何帮助。

英文:

I would like to draw random numbers from a modified exponential distribution:

p(x) = C * a * Exp[-(a*x)^b] with C=1/Gamma[1 + 1/b] for normalization.

How can I do this in julia? Unfortunately I have only little experience with Julia and no experiences with creating custom random numbers. I would be very grateful for any help.

答案1

得分: 6

如果我没有理解错的话,这是一个p-广义高斯分布,在Distributions.jl中有一个相当高效的实现:

using Distributions
mu = 0       # 你的位置参数
alpha = 1/a  # 你的尺度参数
beta = b     # 你的形状参数
p = PGeneralizedGaussian(mu, alpha, beta)

使用Distributions.jl API用于一维分布,你可以通过将其传递给导出的rand方法来从此分布中抽样。以下是从具有mu = 0alpha = 1/2beta = 3PGeneralizedGaussian分布中抽样五个独立标量的示例:

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
  # ... 你的分布参数在这里
end

rand(::AbstractRNG, d::MyDistribution) # 从d中抽样一个值
sampler(d::MyDistribution)             # 返回d或能更高效地从d中抽样的东西
logpdf(d::MyDistribution, x::Real)     # 计算x处d的概率密度函数的对数
cdf(d::MyDistribution, x::Real)        # 计算x处d的累积分布函数
quantile(d::MyDistribution, q::Real)   # 计算d的q分位数
minimum(d::MyDistribution)             # 返回d支持的最小值
maximum(d::MyDistribution)             # 返回d支持的最大值
insupport(d::MyDistribution, x::Real)  # 查询x是否被d支持

该包的文档非常好,所以如果你想学习Julia,这是一个很好的入门方式。

英文:

If I'm not mistaken, that is a p-Generalized Gaussian distribution, which has a rather efficient implementation in 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)

Using the Distributions.jl API for univariate distributions, you can sample from this distribution by passing it to the exported rand method. Here is an example of how to sample five independent scalars from a PGeneralizedGaussian distribution with mu = 0, alpha = 1/2 and 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

If you want to try to implement the distribution yourself, which I do not recommend unless you are doing this as an exercise in programming math in Julia, you need to define a type that holds the static parameters of the distribution (in your case, just the shape parameter and the scale parameter), then define and export the methods listed here to extend the Distributions.jl API to your custom distribution. In particular, you need to define:

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

The documentation of the package is very good, so it's an excellent way to get your feet wet if you are trying to learn Julia.

huangapple
  • 本文由 发表于 2023年6月2日 02:21:23
  • 转载请务必保留本文链接:https://go.coder-hub.com/76384683.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定