模拟均值 = 1 和方差 =0.4 的左截断 Gamma 分布

Question

假设X~Γ(α, β)，我想截断X的所有值

MATLAB 代码：

t = 0.5;  theta = 0.4; 
syms alpha beta 
EX = beta*( igamma(alpha+1,t/beta) /  igamma(alpha,t/beta) ); %Mean  
EX2 = beta^2*( igamma(alpha+2,t/beta) /  igamma(alpha,t/beta) );%Second moment         
VarX = EX2 -EX^2; %Variance       
cond1 = alpha > 0; cond2 = beta > 0; cond3 = EX==1; cond4 = VarX ==theta;
conds =[cond1 cond2 cond3 cond4]; vars = [alpha, beta];
sol = solve(conds, [alpha beta], 'ReturnConditions',true);
soln_alpha = vpa(sol.alpha)
soln_beta = vpa(sol.beta)

以上代码returns仅在放宽α>0的约束条件下才为数字答案。数字答案的 α 为负值，这是错误的，因为 α（形状参数）和 β（比例参数）都必须严格为正。

Answer 1

根据你的标题，我认为你想从均值 = 1 和方差 = 0.4 的 Gamma 分布 生成样本，但希望将分布截断为 [0, inf].

如果X ~ Gamma(alpha,beta), 那么根据定义它必须非负（参见 Gamma 分布 wiki, or MATLAB page). Indeed, both shape and scale parameters are also nonnegative. Note: MATLAB uses the (k,theta) parameterization found on the wiki page.

MATLAB 实现了 probability distribution objects，从从业者（或任何使用数值方法的人）的角度来看，这使得很多事情变得非常方便。

alpha = 0.4;
beta = 0.5;
pd = makedist('Gamma',alpha,beta)   % Define the distribution object

生成样本现在非常容易。

n = 1000;                           % Number of samples
X = random(pd,n,1);                 % Random samples of X ~ Gamma(alpha,beta)

剩下的就是确定形状和尺度参数，使得 E[X] = 1 和 Var(X ) = 0.4。

你需要解决

alpha * beta = E[X],
alpha * (beta^2) = Var(X),

对于 alpha 和 beta。它是具有两个未知数的两个非线性方程组。

但是，截断使这些过时但数值方法可以正常工作。

LB = 0.5;                           % lower bound    (X > LB)
UB = inf;                           % upper bound    (X < UB)
pdt = truncate(pd,LB,UB)            % Define truncated distribution object
Xt = random(pd,n,1);

pdt =
GammaDistribution

Gamma distribution
a = 0.4
b = 0.5
Truncated to the interval [0.5, Inf]

幸运的是，分布 object 的均值和方差无论是否被截断都是可访问的。

mean(pdt)     % compare to mean(pd)
var(pdt)      % compare to var(pd)

您可以用数值方法解决 这个问题来获得您的参数，例如fmincon。

tgtmean = 1;
tgtvar = 0.4;
fh_mean =@(p) mean(truncate(makedist('Gamma',p(1),p(2)),LB,UB));
fh_var =@(p) var(truncate(makedist('Gamma',p(1),p(2)),LB,UB));
fh =@(p) (fh_mean(p)-tgtmean).^2 + (fh_var(p)-tgtvar).^2;
[p, fval] = fmincon(fh,[alpha;beta],[],[],[],[],0,inf)

您可以测试验证的答案：

pd_test = truncate(makedist('Gamma',p(1),p(2)),LB,UB);
mean(pd_test)
var(pd_test)

ans = 1.0377
ans = 0.3758

请注意，由于所需的截断和目标均值，这似乎 ill-conditioned。根据您的应用程序，这可能就足够了。

histogram(random(pd_test,n,1))   % Visually inspect distribution

均值和方差组合必须在基本分布（此处为 Gamma 分布）下可行，但如果截断，则会进一步限制可行参数集。例如，不可能将 X~Gamma() 截断到区间 [5, 500] 并寻求获得 2 的平均值或 600 的平均值。

使用版本 R2017a 验证的 MATLAB 代码。

另请注意，fmincon 等非线性求解器的解决方案可能对某些问题的初始起点敏感。如果该数值方法出现问题，则可能是可行性问题（如上所述），或者可能需要使用多个起点和多个 fmincon 调用来获得多个答案，然后使用最佳答案。

模拟均值 = 1 和方差 =0.4 的左截断 Gamma 分布

Simulate left Truncated Gamma distribution with mean = 1 and variance =0.4

random

matlab

gamma-distribution