如何在Julia中从自定义修正指数分布生成随机数？

对于你提到的修正指数分布（PDF为 ( p(x) = C \cdot a \cdot \exp\left(-(a x)^b\right) )，其中 ( C = 1/\Gamma(1+1/b) )），最适合新手的方法是逆变换采样法，不需要复杂的自定义RNG框架，直接基于均匀分布随机数推导即可。

核心原理：逆变换采样 #

首先推导该分布的累积分布函数（CDF）：

对PDF从0到x积分后，可得 ( F(x) = 1 - \exp\left(-(a x)^b\right) )（归一化常数C的作用刚好让积分结果简化成这个形式）

逆变换采样的逻辑很简单：

1. 生成一个[0,1)区间的均匀随机数 ( u )
2. 解CDF的逆方程得到x：
   ( u = 1 - \exp\left(-(a x)^b\right) )
   整理后得到：( x = \frac{(-\ln(u))^{1/b}}{a} )
   （注：因为1-u也是[0,1)均匀分布，所以直接用u代替1-u不影响结果，计算更简便）

Julia实现代码 #

1. 导入基础包 #

只需要Random生成随机数，SpecialFunctions可选（用于验证PDF的Gamma函数）：

using Random, SpecialFunctions

2. 单个随机数生成函数 #

function modified_exponential_rand(a::Real, b::Real; rng::AbstractRNG=Random.default_rng())
    u = rand(rng)  # 生成[0,1)均匀随机数
    return (-log(u))^(1/b) / a
end

3. 批量生成随机数 #

如果需要生成n个样本，用广播版本效率更高：

modified_exponential_rand(n::Int, a::Real, b::Real; rng::AbstractRNG=Random.default_rng()) = 
    (-log.(rand(rng, n))).^(1/b) ./ a

4. 验证样本正确性 #

可以通过绘制样本直方图和理论PDF对比，确认生成的样本符合分布：

using Plots

# 设置参数
a = 2.0
b = 1.5
sample_count = 100000

# 生成样本
samples = modified_exponential_rand(sample_count, a, b)

# 计算理论PDF
C = 1 / gamma(1 + 1/b)
theoretical_pdf(x) = C * a * exp(-(a*x)^b)

# 绘图对比
histogram(samples, normalize=:pdf, bins=50, label="样本直方图", alpha=0.7)
plot!(theoretical_pdf, 0, maximum(samples), linewidth=2, label="理论PDF", color=:red)
xlabel!("x")
ylabel!("概率密度")
title!("修正指数分布样本验证")

注意事项 #

• 必须保证参数a>0、b>0，否则分布无定义
• 如果需要固定随机种子，在生成样本前调用Random.seed!(your_seed)（比如Random.seed!(1234)），确保结果可复现
• 批量生成时优先用广播版本，比循环效率高很多

内容的提问来源于stack exchange，提问作者HighwayJohn