Расщепление смеси

[Александрович]

Emergency в сообщении #1436474 писал(а):

я вижу только 19 точек.

Ещё раз.

Александрович в сообщении #1436344 писал(а):

Emergency в сообщении #1436340 писал(а):

И точек всего 19, а не 196.

Это не точек, а интервалов. Например модальный интервал $[89;112]$ содержит 51 точку.

-- Чт янв 23, 2020 04:37:47 --

Евгений Машеров в сообщении #1436475 писал(а):

Как-то очень уж похоже на логнормальное...

Lia в сообщении #1436477 писал(а):

Или хи-квадрат (?)

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

[Emergency]

Александрович в сообщении #1436483 писал(а):

Ещё раз.

Спасибо, я уже разобрался, что график - это не результат эксперимента, а его статобработка.

[Александрович]

Emergency в сообщении #1436527 писал(а):

график - это не результат эксперимента

На графике серые столбики это результаты реальных наблюдений.

[Евгений Машеров]

Ну, вообще-то оценивание минимизацией хи-квадрат метод известный. Кендалл и Стьюарт, "Статистические выводы и связи", 19.25 - 19.28 (стр. 132-136)
Или Вы именно применительно к разделению смесей?

[Александрович]

Евгений Машеров в сообщении #1436544 писал(а):

оценивание минимизацией хи-квадрат метод известный

А правомерно ли после нахождения параметров этим методом проверять гипотезу о принадлежности выборочного распределения гипотетическому по критерию хи-квадрат?

[B@R5uk]

Александрович в сообщении #1436197 писал(а):

Собрал реальные 196 значений некой случайной величины.

А вы можете привести сырые результаты в виде таблицы чисел поиграться?

[Александрович]

Используется синтаксис Text

29,4
40,5
42,5
42,7
45,3
45,8
47,2
47,4
48,2
50,9
50,9
56,5
56,8
58,3
59,1
59,1
59,9
60,4
60,6
60,8
61,1
61,4
62,7
62,8
63,3
64,1
64,3
64,7
66,1
66,5
67,0
67,0
67,5
67,9
68,2
68,3
69,3
70,6
71,7
71,9
72,2
73,3
73,3
73,9
74,8
74,8
75,2
76,5
76,9
77,0
77,2
77,3
77,3
77,4
77,4
77,7
78,9
78,9
79,4
79,6
79,7
79,8
80,5
81,0
81,3
82,0
82,0
82,4
83,1
83,3
83,4
83,9
84,9
84,9
86,0
87,0
87,1
88,5
89,1
89,6
89,6
89,6
89,8
90,0
90,1
90,8
90,9
91,8
92,9
93,0
93,1
93,2
93,7
94,0
94,0
94,2
95,2
95,7
95,9
96,0
96,0
96,9
97,9
99,5
99,7
100,1
100,3
100,5
100,7
104,3
104,4
104,8
105,1
105,3
105,3
105,6
105,8
106,1
106,8
107,3
107,7
108,6
108,8
109,8
109,8
110,3
110,3
110,9
111,0
112,2
113,7
113,9
115,5
115,5
116,0
116,4
117,3
117,5
117,8
119,0
119,3
119,3
120,0
120,5
121,0
121,6
124,0
124,7
124,9
127,3
127,9
130,9
131,7
133,9
136,5
136,8
143,1
145,2
145,8
147,5
148,8
150,5
151,3
155,2
156,1
158,0
159,1
160,8
161,3
164,7
166,8
171,6
172,6
173,8
175,4
190,3
191,1
196,5
200,4
200,7
203,3
203,6
204,8
206,2
209,0
209,3
215,7
223,4
224,2
234,9
240,5
258,5
285,8
325,1
403,3
454,1

[B@R5uk]

Александрович, спасибо. И спрячьте под спойлер, пока не поздно. Спасибо.

За этим распределением есть какая-либо физическая теория, что вы выбрали сумму двух нормальных? Есть ли какая-либо объективная причина, что нет значений меньшей величины (например, щель блокировала вылет частиц и значения с меньшими величинами механически отбрасывались)?

[Александрович]

B@R5uk в сообщении #1436583 писал(а):

выбрали сумму двух нормальных

Это значение размаха виброперемещений за 196 оборотов вращения ротора. Полагаю что источников порождения колебаний по меньшей мере два. А нормальные потому что можно представить в виде суммы многих случайных факторов.

[B@R5uk]

Евгений Машеров в сообщении #1436475 писал(а):

Как-то очень уж похоже на логнормальное...

Это если раскидать случайные значения по бинам и построить дискретную плотность вероятности. Если же посмотреть на честную функцию распределения, то не очень то и похоже:


Кстати, в логарифмическом масштабе выглядит красивше:


-- 23.01.2020, 18:45 --

Александрович в сообщении #1436587 писал(а):

Полагаю что источников порождения колебаний по меньшей мере два. А нормальные потому что можно представить в виде суммы многих случайных факторов.

Вообще, тот факт, что распределение X является суммой двух других распределений A и B с весами w и 1-w, означает что события из распределений A и B происходят исключающе. То есть либо одно, либо другое. Это слишком жёсткая взаимозависимость для двух источников случайных колебаний. Но я покритиковать покритиковал, а предложить взамен чего-то более разумного не могу. Потому что не знаю. Надеюсь, более опытные коллеги подскажут.

За то, что я знаю наверняка — это то, что искать параметры выбранного распределения по заданной выборке надо через функцию правдоподобия. Разбиение на бины — это искусственный приём и приводит к большой потере информации.

[B@R5uk]

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



Хвост распределения жутко плохо зафитился. Надо бы что-нибудь другое взять.

-- 23.01.2020, 20:11 --

Хотя нет. Если выкинуть из расчёта последние четыре точки, то фитится замечательно:

Видимо у гаусса хвост слишком быстро спадает, по сравнению с экспериментом, поэтому точки в хвосте дают большие искажения для данного фита.

Наверное, стоит ещё попробовать суммой двух гамма-распределений зафитить. В вики вон чё пишут:

Occurrence and applications писал(а):

The gamma distribution has been used to model the size of insurance claims[16] and rainfalls.[17] This means that aggregate insurance claims and the amount of rainfall accumulated in a reservoir are modelled by a gamma process – much like the exponential distribution generates a Poisson process.

[EUgeneUS]

Александрович в сообщении #1436587 писал(а):

Это значение размаха виброперемещений за 196 оборотов вращения ротора.

Все таки нужно больше информации, как это измеряется, и почему не бывает отрицательных значений.
Ниже предположения и домыслы для примера.
1. Что-то вращается и из-за вибрации как-то смещается.
2. Перемещение измеряется относительно чего? Вероятно, относительно какой-то станины, на которой стоит механизм, а значит перемещение происходит по двум осям.
3. Есть две существенно разные оси - вдоль и поперек оси вращения.
4. Пусть проекция смещения по каждой оси имеет центральное нормальное распределение, но с различной дисперсией.
5. Тогда распределение модуля перемещения будет нецентральным хи-распределением (без квадрата) с двумя степенями свободы.

[B@R5uk]

Попробовал зафитить суммой двух гамма-распределений. На мой взгляд результат стал значительно лучше. Это может на графике не особо заметно, но в том месте, где густота точек значительна теоретическая функция распределения значительно лучше повторяет форму практических значений. Правда, в формулах я использовать не стандартную параметризацию, а такую (чтобы лучше соответствовало параметризации гауссианины):

$$f\left( x,\mu ,\sigma \right)={{\sigma }^{-1}}{{\Gamma }^{-1}}\left( \frac{\mu }{\sigma } \right)\exp \left( -\frac{x}{\sigma }+\left( \frac{\mu }{\sigma }-1 \right)\ln \frac{x}{\sigma } \right)$$

$$F\left( x,k,\sigma \right)={{\Gamma }^{-1}}\left( \frac{\mu }{\sigma } \right)\gamma \left( \frac{\mu }{\sigma },\frac{x}{\sigma } \right)$$

Результат работы такой:



Если кто захочет с программой поиграться, то вот код:

Используется синтаксис Matlab M

clc
clearvars
format compact

%   https://dxdy. ru/topic138642.html

%fit = 'gauss';
fit = 'gamma';

wdata = [
    29.4    40.5    42.5    42.7    45.3    45.8    47.2
    47.4    48.2    50.9    50.9    56.5    56.8    58.3
    59.1    59.1    59.9    60.4    60.6    60.8    61.1
    61.4    62.7    62.8    63.3    64.1    64.3    64.7
    66.1    66.5    67.0    67.0    67.5    67.9    68.2
    68.3    69.3    70.6    71.7    71.9    72.2    73.3
    73.3    73.9    74.8    74.8    75.2    76.5    76.9
    77.0    77.2    77.3    77.3    77.4    77.4    77.7
    78.9    78.9    79.4    79.6    79.7    79.8    80.5
    81.0    81.3    82.0    82.0    82.4    83.1    83.3
    83.4    83.9    84.9    84.9    86.0    87.0    87.1
    88.5    89.1    89.6    89.6    89.6    89.8    90.0
    90.1    90.8    90.9    91.8    92.9    93.0    93.1
    93.2    93.7    94.0    94.0    94.2    95.2    95.7
    95.9    96.0    96.0    96.9    97.9    99.5    99.7
   100.1   100.3   100.5   100.7   104.3   104.4   104.8
   105.1   105.3   105.3   105.6   105.8   106.1   106.8
   107.3   107.7   108.6   108.8   109.8   109.8   110.3
   110.3   110.9   111.0   112.2   113.7   113.9   115.5
   115.5   116.0   116.4   117.3   117.5   117.8   119.0
   119.3   119.3   120.0   120.5   121.0   121.6   124.0
   124.7   124.9   127.3   127.9   130.9   131.7   133.9
   136.5   136.8   143.1   145.2   145.8   147.5   148.8
   150.5   151.3   155.2   156.1   158.0   159.1   160.8
   161.3   164.7   166.8   171.6   172.6   173.8   175.4
   190.3   191.1   196.5   200.4   200.7   203.3   203.6
   204.8   206.2   209.0   209.3   215.7   223.4   224.2
   234.9   240.5   258.5   285.8   325.1   403.3   454.1
];

switch fit
    case 'gauss'
        thrown_num = 4;
        rsqrt2pi = 1 / sqrt(2 * pi);
        base_density = @(x, m, s)...
            rsqrt2pi * exp(- (x - m) .^ 2 / (2 * s ^ 2)) / s;
        base_distrib = @(x, m, s)...
            (1 + erf((x - m) / (sqrt(2) * s))) / 2;
        title_str = 'Фит экспериментальных данных суммой двух нормальных распределений';
    case 'gamma'
        thrown_num = 3;
        base_density = @(x, m, s)...
            exp(- x / s + (m / s - 1) * log(x / s)) / (s * gamma(m / s));
        base_distrib = @(x, m, s) gammainc(x / s, m / s);
        title_str = 'Фит экспериментальных данных суммой двух гамма-распределений';
end

wdata = wdata';
wdata = wdata(1 : end - thrown_num);
dsize = numel(wdata);

final_density = @(x, w, m1, s1, m2, s2)...
    w * base_density(x, m1, s1) + (1 - w) * base_density(x, m2, s2);
final_distrib = @(x, w, m1, s1, m2, s2)...
    w * base_distrib(x, m1, s1) + (1 - w) * base_distrib(x, m2, s2);

taget_func = @(z) -sum(log(final_density(wdata,...
    min(1, max(0, z(1))), z(2), z(3), z(4), z(5))));

z0 = [0.6, 90, 8, 170, 12];

fmso = optimset('fminsearch');
fmso = optimset(fmso, 'Display', 'final', 'TolFun', 1e-12, 'TolX', 1e-12,...
    'MaxIter', 20000, 'MaxFunEvals', 10000);
z = fminsearch(@(z) taget_func(z), z0, fmso);

%z(5) = 55;
%z = z0;

%disp(final_density(100, z0(1), z0(2), z0(3), z0(4), z0(5)))
%disp(final_distrib(500, z0(1), z0(2), z0(3), z0(4), z0(5)))

xx = 0 : 0.1 : 300;
yy = final_distrib(xx, z(1), z(2), z(3), z(4), z(5));

subplot(1, 1, 1)
plot(wdata, (1 : dsize) / dsize, '>');
hold on
plot(xx, yy, 'k', 'LineWidth', 2)
hold off
xlim([0, 300])
ylabelvalues = 0 : 0.2 : 1;
set(gca, 'YTick', ylabelvalues)
set(gca, 'YTickLabel', ylabelvalues)
grid on
title(title_str)
xlabel('Случайная величина')
ylabel('Функция распределения')
legend('Эксперимент', 'Теория', 'Location', 'NorthWest')
text(240, 0.355, {['w = ', num2str(z(1))],...
    ['\mu_1 = ', num2str(z(2))], ['\sigma_1 = ', num2str(z(3))],...
    ['\mu_2 = ', num2str(z(4))], ['\sigma_2 = ', num2str(z(5))]})

%semilogx(wdata, 1 : dsize, '.');
%xlim([20, 500])
%xlabelvalues = [20, 30, 40, 60, 80, 100, 150, 200, 300, 400, 500];
%set(gca, 'XTick', xlabelvalues)
%set(gca, 'XTickLabel', xlabelvalues)
 

[Emergency]

EUgeneUS в сообщении #1436616 писал(а):

Перемещение измеряется относительно чего? Вероятно, относительно какой-то станины, на которой стоит механизм

А колебания станины со своей резонансной частотой не могут дать второй источник смещений?

[Александрович]

Emergency, вполне возможно.

-- Пт янв 24, 2020 03:42:43 --

EUgeneUS в сообщении #1436616 писал(а):

почему не бывает отрицательных значений

Размах виброперемещения в принципе не может быть отрицательным.