Имитация левого усеченного гамма-распределения со средним значением = 1 и дисперсией = 0,4 - PullRequest
Новая
       46

Имитация левого усеченного гамма-распределения со средним значением = 1 и дисперсией = 0,4

0 голосов
/ 26 апреля 2019

Предположим, X ~ & Gamma; (& alpha ;, & beta;), я хотел бы усечь все значения 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)

Приведенный выше код возвращает числовой ответ, только если ограничение, что & alpha;> 0 ослаблено. Числовой ответ имеет отрицательное значение & alpha; что неверно, так как оба & alpha; (параметр формы) и & beta; (параметр масштаба) должен быть строго положительным.

1 Ответ

1 голос
/ 26 апреля 2019

Исходя из вашего названия, я предполагаю, что вы хотите сэмплировать из гамма-распределения со средним значением = 1 и дисперсией = 0,4, но хотите, чтобы распределение было усечено до [0, inf] .

Если X ~ Gamma ( alpha , beta ), то по определению он должен быть неотрицательным (см. Распределение гаммы wiki , или страница MATLAB ). Действительно, параметры формы и масштаба также неотрицательны. Примечание: MATLAB использует ( k , theta ) параметризацию, найденную на вики-странице .

В MATLAB реализовано объектов распределения вероятностей , которые делают многие вещи очень удобными с точки зрения практиков (или тех, кто использует численные подходы). 

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.

Вам нужно решить

альфа * бета = E [ X ],
alpha * ( beta ^ 2) = Var ( X ),

для альфа и бета . Это система из двух нелинейных уравнений с двумя неизвестными.

Однако, усечение делает эти устаревшие, но численные подходы, будут работать нормально. 

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 =
Гаммараспределение

Гамма-распределение
а = 0,4
б = 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

Histogram of Gamma dist truncated to [0.5, inf]

Комбинации среднего значения и дисперсии должны быть возможными при базовом распределении (здесь, гамма-распределение), но в случае усечения это еще больше ограничивает набор возможных параметров. Например, было бы невозможно усечь X ~ Gamma () до интервала [5, 500] и попытаться получить среднее значение 2 или среднее значение 600.

Код MATLAB подтвержден версией R2017a.

Также обратите внимание, что решения от нелинейных решателей, таких как fmincon, могут быть чувствительными к начальной точке для некоторых задач. Если при таком численном подходе возникают проблемы, это может быть проблемой осуществимости (как упоминалось выше) или может потребоваться использование нескольких начальных точек и нескольких вызовов fmincon для получения нескольких ответов, а затем используйте лучший.

...