使用mean = 1和variance = 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)

仅当放宽α> 0的约束时，上述代码才返回数字答案。数字答案具有α的负值，这是错误的，因为α（形状参数）和β（比例参数）必须严格为正。

Answer 1

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

如果X~Gamma（alpha，beta），则根据定义它必须是非负的（参见Gamma Distribution wiki，或MATLAB page）。实际上，形状和尺度参数也都是非负的。注意：MATLAB使用wiki page上的（k，theta）参数化。

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 = a [k]， Alpha *（beta ^ 2）=次（ks），

对于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 = 伽玛分布

Gamma分布 a = 0.4 b = 0.5 截断到区间[0.5，Inf]

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

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

请注意，由于所需的截断和目标均值，这似乎是病态的。根据您的应用，这可能已经足够了。

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

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

---

使用R2017a版验证的MATLAB代码。

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