Как исправить размер выходной выборки выборки асимметричного трапециевидного распределения

Я пытаюсь создать перекошенное трапециевидное распределение, используя выборку обратного преобразования. Входными данными являются значения, где начинаются и заканчиваются линейные изменения (a, b, c, d), а также размер выборки.

a=-3;b=-1;c=1;d=8; 
SampleSize=10e4;
h=2/(d+c-a-b);

Затем я рассчитываю соотношение длины пандусов и плоских компонентов, чтобы получить размер выборки для каждого:

firstramp=round(((b-a)/(d-a)),3);
flat=round((c-b)/(d-a),3);
secondramp=round((d-c)/(d-a),3);
n1=firstramp*SampleSize; %sample size for first ramp
n3=secondramp*SampleSize; %sample size for second ramp
n2=flat*SampleSize;

И затем, наконец, я получаю гистограмму из следующего кода:

quartile1=h/2*(b-a);
quartile2=1-h/2*(d-c);

y1=linspace(0,quartile1,n1);
y2=linspace(quartile1,quartile2,n2);
y3=linspace(quartile2,1,n3);

%inverse cumulative distribution functions
invcdf1=a+sqrt(2*(b-a)/h)*sqrt(y1);
invcdf2=(a+b)/2+y2/h;
invcdf3=d-sqrt(2*(d-c)/h)*sqrt(1-y3);

distr=[invcdf1 invcdf2 invcdf3];

histogram(distr,100)

Однако выборка пандусов и плоских компонентов неодинакова, выглядит так:

Неверный вывод

Я исправил это методом проб и ошибок, уменьшив размер выборки рамп вдвое:

n1=0.5*firstramp*SampleSize; %sample size for first ramp
n3=0.5*secondramp*SampleSize; %sample size for second ramp
n2=flat*SampleSize;

В результате распределение выглядело так:

Фиксированное распространение

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

Я также пробовал разные комбинации изменения размеров сэмплов рамп и флэт. Это также работает:

n1=0.75*firstramp*SampleSize; %sample size for first ramp
n3=0.75*secondramp*SampleSize; %sample size for second ramp
n2=1.5*flat*SampleSize;

Это увеличивает выходные образцы, но это все еще не близко.

Любая помощь будет оценена.

Полный код:

a=-3;b=-1;c=1;d=8; 
SampleSize=10e4;%*1.33333333333333;
h=2/(d+c-a-b);
firstramp=round(((b-a)/(d-a)),3);
flat=round((c-b)/(d-a),3);
secondramp=round((d-c)/(d-a),3);

n1=firstramp*SampleSize; %sample size for first ramp
n3=secondramp*SampleSize; %sample size for second ramp
n2=flat*SampleSize;

quartile1=h/2*(b-a);
quartile2=1-h/2*(d-c);

y1=linspace(0,quartile1,.75*n1);
y2=linspace(quartile1,quartile2,1.5*n2);
y3=linspace(quartile2,1,.75*n3);

%inverse cumulative distribution functions
invcdf1=a+sqrt(2*(b-a)/h)*sqrt(y1);
invcdf2=(a+b)/2+y2/h;
invcdf3=d-sqrt(2*(d-c)/h)*sqrt(1-y3);

distr=[invcdf1 invcdf2 invcdf3];

histogram(distr,100)
%end

matlab algorithm math distribution skew
person Sumit Mann    schedule 16.12.2020    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Я не знаю Matlab, поэтому я надеялся, что кто-то еще включится в это, но, поскольку здесь никто этого не сделал.

Если я правильно читаю ваш код, то, что вы сделали, не является инверсией. Инверсия 1-1, т. е. один единый вход дает один результат. Вы, кажется, используете технику, известную как метод композиции. По составу общий дистрибутив состоит из составных частей, каждую из которых легко сгенерировать. Вы выбираете, из какого компонента генерировать, исходя из их пропорций/вероятностей по отношению к целому. Для функций плотности вероятность определяется как площадь под кривой плотности, поэтому ваша первая ошибка заключалась в выборке компонентов относительно ширины каждого компонента, а не в использовании их площадей. Правильные пропорции выборки: 2/13, 4/13 и 7/13 для компонентов firstramp, flat и secondramp соответственно. Вторая ошибка (относительно незначительная) заключалась в назначении точных размеров выборки для каждого из компонентов. Наличие вероятности 2/13 не означает, что ровно 2*SampleSize/13 ваших выборок будут из firstramp, это означает, что это ожидаемый размер выборки для этого компонента. Ожидаемое значение случайной переменной не обязательно (и даже не может быть) тем результатом, который вы действительно получите.

В псевдокоде композиционный подход будет

generate U ~ Uniform(0,1)
if U <= 2/13:
   generate and return a value from firstramp
else if U <= 6/13:
   generate and return a value from flat
else:
   generate and return a value from secondramp

Обратите внимание, что поскольку каждый из вариантов generate будет использовать один или несколько униформ, а выбор между вариантами требует униформы U, это не инверсия.

Если вам нужна реальная инверсия, вам нужно количественно определить свою плотность, проинтегрировать ее, чтобы получить кумулятивную функцию распределения, а затем применить метод инверсии, установив F(X) = U и найдя X. Поскольку ваше распределение состоит из отдельных компонентов, плотность и кумулятивная плотность будут кусочными функциями.

После получения высоты на основе требования о том, что площади двух треугольников и плоского участка должны составлять в сумме 1, я пришел к следующему для вашей плотности:

       | (x + 3) / 13       -3 <= x <= -1
       |
f(x) = | 2 / 13             -1 <= x <= 1
       |
       | 2 * (8 - x) / 91    1 <= x <= 8

Интеграция этого и сбор терминов дает CDF:

       | (x + 3)**2 / 26                    -3 <= x <= -1
       |
F(x) = | (2 + x) * 2 / 13                   -1 <= x <= 1
       |
       | 6 / 13 + [49 - (x - 8)**2] / 91     1 <= x <= 8

Наконец, определение значений F(x) в точках разрыва между сегментами и применение инверсии дает следующий алгоритм псевдокода:

generate U ~ Uniform(0,1)
if U <= 2 / 13:
    return 2 * sqrt( (13 * U) / 2 ) - 3
else if U <= 6 / 13:
    return (13 * U) / 2 - 2:
else:
    return 8 - sqrt( 91 * (1 - U) )

Обратите внимание, что это истинная инверсия. Результат определяется созданием одного U и преобразованием его различными способами в зависимости от того, в какой диапазон он попадает.

person pjs    schedule 30.12.2020