Критерий согласия для распределения Парето

В математической статистике значимость исследования выборки очень существенна. Выборка подвергается обработке и выдвигается предположение о распределении, которому подчиняется выборка.

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

1. Теоретическая часть

1.1 Распределение Парето

Распределение Парето — это двухпараметрическое семейство абсолютно непрерывных распределений.

Функция распределения F(x) имеет вид (1.1).

(1.1)

где б, x₀ — параметры распределения, x > x₀ > 0, б > 0.

Функция плотности распределения f(x) имеет вид (1.2).

(1.2)

График функции распределения приведен на рисунке 1. График функции плотности распределения приведен на рисунке 2.

Рисунок 1 — График функции распределения Рисунок 1 — График функции плотности распределения

1.2 Критерий согласия ч²

Предположим, что по виду гистограммы или полигона частот или из каких-либо других соображений удается выдвинуть гипотезу о множестве функций определенного вида (нормальных, показательных, биномиальных и т. п.), к которому может принадлежать функция распределения исследуемой случайной величины X. Критерий ч² Пирсона (критерий согласия ч²) позволяет производить проверку согласия эмпирической функции распределения F*(x) с гипотетической функцией распределения F (x).

Для этого придерживаются следующей последовательности действий:

1) Диапазон изменения экспериментальных данных разбивается на k интервалов;

2) На основании гипотетической функции F(x) вычисляют вероятность попадания с.в. X в частичные интервалы [x_i-1, x_i] по формуле (1.3);

p_i = P (x_i-1 ? X? x_i), i=1,2., k (1.3)

3) Умножая полученные вероятности p_i на объем выборки n, получают теоретические частоты np_i частичных интервалов [x_i-1, x_i], т. е. частоты, которые следует ожидать, если гипотеза справедлива;

4) Вычисляю выборочную статистику (критерий) ч² по формуле (1.4).

(1.4)

где m_i — количество значений с.в., попавших в i-й интервал;

n — объём выборки.

Если гипотеза верна, то при n>? распределение выборочной статистики, независимо от вида функции F (x), стремится к распределению ч² с v= k-r-1 степенями свободы (k — число частичных интервалов, r — число параметров гипотетической функции F(x), оцениваемых по данной выборке).

Критерий ч² сконструирован таким образом, что чем ближе к нулю наблюдаемое значение критерия ч², тем вероятнее, что гипотеза справедлива. Поэтому для проведения гипотезы применяется критерий ч² с правосторонней критической областью. Необходимо найти по таблице квантилей ч²-распределения по заданному уровню значимости б и числу степеней свободы v критическое значение ч²_б,v, удовлетворяющее условию p (ч²? ч²_б,v) = б.

Если ч²_набл.? ч²_б,v, то считается, что гипотетическая функция F (x) не согласуется с результатами эксперимента. Если ч²_набл.? ч²_б,v, то считается, что гипотетическая функция F (x) согласуется с результатами эксперимента.

1.3 Алгоритм обработки выборки

1) Сортируем выборку по возрастанию (преобразуем в вариационный ряд)

2) Находим минимальный x_min и максимальный x_max элемент выборки

3) Находим длину интервалов группировки h по формуле (1.5)

(1.5)

где k — число интервалов группировки.

4) Находим левые x_l и правые x_r границы интервалов группировки по формулам (1.6)

(1.6)

5) Находим центры x_k^* интервалов группировки по формуле (1.7)

(1.7)

6) Для каждого интервала группировки (x_k-1, x_k) находим число n_k^* (абсолютная частота) элементов выборки, попавших в этот интервал. Важно чтобы каждый элемент выборки был отнесен к одному и только к одному интервалу, а если значение элемента попадает на границу интервала, то его относят к интервалу с младшим номером. Минимальный элемент всегда относится к первому интервалу, максимальный к последнему.

7) Вычисляем относительные частоты Otn_k^* по формуле (1.8) как отношение абсолютной частоты к объему выборки. Убеждаемся, что сумма всех относительных частот равна единице (допускается небольшое отличие от единицы в рамках погрешности вычислений).

(1.8)

8) Строим гистограмму относительных частот — фигуру, состоящую из k прямоугольников, опирающихся на интервалы группировки. Площадь k-го прямоугольника полагают равной относительной частоте данного интервала. Высота k-го прямоугольника H_k рассчитывается по формуле (1.9).

(1.9)

Убеждаемся, что сумма всех высот H_k, умноженная на h, равна единице (допускается небольшое отличие от единицы в рамках погрешности вычислений). На оси абсцисс выбираем начальную точку чуть левее точки x_min, и такой масштаб, чтобы на оси поместился интервал [x_min, x_max] и отчетливо различались точки x_l, x_r. На оси ординат выбираем начало отсчета в точке 0 и такой масштаб, чтобы отчетливо различались H_k. Для построения гистограммы относительных частот на ось абсцисс наносим интервалы [x_l, x_r] и, используя каждый из них как основание, строим прямоугольник с соответствующей высотой H_k. Получаем гистограмму.

9) Вычисляем параметры распределения Парето б и x₀. Для этого используем систему (1.10).

(1.10)

где, ,, .

Получается система (1.11):

(1.11)

Решая систему уравнений, находим параметры б и x₀. При этом должны выполняться условия б > 0. Параметры распределения вычисляются по формулам (1.12).

(1.12)

10) Строим график функции плотности распределения f (x) с вычисленными параметрами распределения б и x₀, где x — это значения центров интервалов x_k*.

11) Проверяем выборку на соответствие распределению Парето по критерию согласия ч² по формуле (1.4), где m_i = n_k^*, n - объем выборки, p_i = F(x_r) — F(x_l), или по формуле интеграла (p_i = x₀^б((x_l)^-б — (x_r)^-б). При этом минимальное значение по левой границе равно значению больше нуля, а максимальное значение по правой границе — бесконечность ?, а также должны выполняться условия функции распределения. Хи-квадрат крит. ч²_{б,v =} ч²_крит. находится по таблице, где v = 3, б — выбирается из таблицы.

2. Практическая часть

Для автоматизации обработки выборки была разработана программа Pareto_distribution.exe, в которую заложены алгоритмы обработки выборки и возможность быстрого получения результата. Вид программы представлен на рисунке 3.

парето распределение программа выборка Рисунок 3 — Программа Pareto_distribution.exe

В программе заложена возможность ввода либо неупорядоченной выборки, либо статистического ряда. Для неупорядоченной выборки возможен импорт из внешнего файла с расширением *.xml и *.txt. Также пользователь может сам вводить выборку.

Для разбиения выборки на интервалы заложено определенное число интервалов, на которые можно разбить выборку. Это число интервалов соответствует числу интервалов из таблицы ч²_крит.

Для контроля правильности обработки выборки и лучшего понимания самого процесса обработки представлена таблица, в которой отображаются все вычисляемые данные, заложенные в постановке задачи.

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

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

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

Для импорта выборки из файла существуют некоторые правила:

1) Для текстовых файлов:

— выборка в поле записывается в один столбик, начиная каждое значение с новой строки;

— в импортируемом файле выборка должна быть записана в один столбик, каждое значение с новой строки. В файле не должно присутствовать лишних значений и заголовков. Пример файла с выборкой представлен на рисунке 4.

Рисунок 4 — Импорт текстового файла

2) Для файлов формата *.xls (MS Excel):

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

— для статистического ряда: выборка представляет собой два столбика, в первый из которых записано значение выборки, а во второй записана частота этого значения. Предполагается, что первая строка — заголовок столбца. Пример статистического ряда представлен на рисунке 5.

Рисунок 5 — Импорт из Excel

При импорте или вводе выборке необходимо, чтобы в таблице в программе не было пустых или не заполненных строк. Если такие строки остаются, необходимо встать на необходимую ячейку и нажать кнопку «Удалить». Пример правильно заполнения таблицы приведен на рисунке 6.

Рисунок 6 — Правильный ввод выборки

2.1 Пример 1

Дана выборка, записанная в файле формата *.txt. Сделать выводы по соответствию данной выборке распределению Парето, используя различные интервалы разбиения и уровень значимости б.

В программе переходим на вкладку «Неупорядоченная выборка (*.txt)» и нажимаем кнопку «Импорт». Появляется окошко, представленное на рисунке 7.

Рисунок 7 — Окно импорта

Выбираем нужный файл и щелкаем кнопку «Open». Выборка, записанная в файле, отобразится в окошечке программы, как показано на рисунке 8.

Рисунок 8 — Выборка из файла После этого с помощью движка выбираем количество интервалов, на которое хотим разбить выборку. Пусть k=10. После этого нажимаем кнопку «Разбить». В таблице отобразятся данные, полученные при разбивке выборки. Полученный результат показан на рисунке 9. С помощью полосы прокрутки можно посмотреть всю таблицу. Для контроля правильности обработки выборки в программу внедрен специальный блок.

Рисунок 9 — Данные разбивки

Для изменения уровня значимости б щелкаем на нужное значение один раз. Выбираем, например первое значение.

Получившиеся гистограмма относительных частот и график функции плотности распределения представлены на рисунке 10.

При данных значениях разбивки делаем вывод, что данная выборка согласуется с распределением Парето, т. е. выполнены условия критерия согласия Пирсона. Вывод указан на рисунке 11.

Для того, чтобы сохранить полученные данные обработки выборки, в программе находится кнопка «Сохранить». При ее нажатии, все полученные данные сохранятся в формате *.xml. Сохраненный файл представлен на рисунке 12.

Рисунок 10 — График Рисунок 11 — Вывод Рисунок 12 — Экспорт в Excel

2.2 Пример 2

Дана выборка в виде статистического ряда. Нужно произвести произвольное разбиение выборки и посмотреть, при какой значении k выборка подчиняется распределению Парето.

В программе переходим на вкладку «Статистический ряд» и нажимаем кнопку «Импорт». Появляется диалоговое окно, с помощью которого открываем файл с готовой выборкой. Нажимаем кнопку «Open».

Рисунок 13 — Статистический ряд

После этого переходим к разбиению выборки:

1) Пусть k=10. Тогда получим данные, приведенные на рисунке 14. По гистограмме относительных частот видно, что выборка не подчиняется распределению Парето. Имеем ч²_крит. > ч²_теор., следовательно, не выполнено условие Пирсона;

Рисунок 14 — Итоги обработки

2) Пусть k=6. Имеем ч²_крит. < ч²_теор., следовательно, условие Пирсона выполнено; а значит, выборка при таком разбиении подчиняется закону распределения Парето.

Заключение

Критерий согласия Пирсона является очень удобным инструментом для определения, относится ли случайная выборка к тому или иному распределению.

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

Приложение А

function Xls_To_StringGrid (AGrid: TStringGrid; AXLSFile: string): Boolean;

const

xlCellTypeLastCell = $ 00B;

var

XLApp, Sheet: OLEVariant;

RangeMatrix: Variant;

SaveChanges: OleVariant;

x, y, k, r: Integer;

i:integer;

begin

Result:= False;

// Create Excel-OLE Object

XLApp:= CreateOleObject ('Excel. Application');

try

// Hide Excel

XLApp. Visible:= False;

// Open the Workbook

XLApp. Workbooks. Open (AXLSFile);

// Sheet:= XLApp. Workbooks[1]. WorkSheets[1];

Sheet:= XLApp. Workbooks [ExtractFileName (AXLSFile)]. WorkSheets[1];

// In order to know the dimension of the WorkSheet, i. e the number of rows

// and the number of columns, we activate the last non-empty cell of it

Sheet. Cells. SpecialCells (xlCellTypeLastCell, EmptyParam).Activate;

// Get the value of the last row

x:= XLApp. ActiveCell. Row;

// Get the value of the last column

y:= XLApp. ActiveCell. Column;

// Set Stringgrid’s row &col dimensions.

AGrid. RowCount:= x;

// Assign the Variant associated with the WorkSheet to the Delphi Variant

RangeMatrix:= XLApp. Range ['A1', XLApp. Cells. Item [X, Y]]. Value;

//XLApp. Quit;

// Define the loop for filling in the TStringGrid

{K:= 1;

repeat

for R:= 1 to Y do

AGrid. Cells[®, (K)]: = RangeMatrix [K+1, R];

AGrid. RowCount:= K + 1;

Inc (K, 1);

until

K >= (X-1);}

K:=1;

While K

for R:= 1 to Y do //begin

AGrid. Cells[®, (K)]: = RangeMatrix [K+1, R];

AGrid. RowCount:= K + 1;

Inc (K, 1);

//end;

end;

For i:=1 to AGrid. RowCount-1 do

AGrid. Cells [0, i]: =IntToStr (i);

// Unassign the Delphi Variant Matrix

RangeMatrix:= Unassigned;

finally

// Quit Excel

if not VarIsEmpty (XLApp) then

begin

XLApp. DisplayAlerts:= False;

XLApp. Quit;

XLAPP:= UnAssigned;

Sheet:= UnAssigned;

Result:= True;

end;

end;

end;

procedure InsertFileInMemo (RzMemo1: TRzMemo; FileName: string;

ReplaceSel: Boolean);

var

Stream: TMemoryStream;

NullTerminator: Char;

begin

Stream:= TMemoryStream. Create;

try

Stream. LoadFromFile (FileName);

Stream. Seek (0, 2);

NullTerminator:= #0;

Stream. Write (NullTerminator, 1);

if not ReplaceSel then

RzMemo1. SelLength:= 0;

SendMessage (RzMemo1. Handle, EM_ReplaceSel, 0,

LongInt (Stream. Memory));

finally

Stream. Free;

end;

end;

procedure Sort;

var i: integer;

changed: boolean;

buf:real;

begin

repeat

Changed:= FALSE;

for i:= 0 to n-2 do

if V[i] > V [i+1] then begin

buf:= V[i];

V[i]:=V [i+1];

V [i+1]: = buf;

changed:= TRUE;

end;

until

not changed;

end;

procedure Vkladka3;

var i: integer;

begin

n:=Form1. RzMemo1. Lines. Count;

SetLength (V, n);

for i:=0 to n-1 do begin

V[i]: =StrToFloat (Form1. RzMemo1. Lines. Strings[i]);

end;

Sort;

end;

procedure Vkladka1;

var i: integer;

begin

n:=Form1.SG. RowCount-1;

SetLength (V, n);

for i:=0 to n-1 do begin

V[i]: =StrToFloat (Form1.SG. Cells [1, i+1]);

end;

Sort;

end;

procedure Vkladka2;

var i, j, m: integer;

l, k, r, w: integer;

SummaEl:integer;

MXi:array of real;

MNi:array of integer;

begin

SummaEl:=0;

m:=Form1.SG1. RowCount-1;

for i:=0 to m-1 do

SummaEl:=SummaEl+StrToInt (Form1.SG1. Cells [2, i+1]);

n:=SummaEl;

SetLength (MXi, m);

SetLength (MNi, m);

for i:=0 to m-1 do begin

MXi[i]: =StrToFloat (Form1.SG1. Cells [1, i+1]);

MNi[i]: =StrToInt (Form1.SG1. Cells [2, i+1]);

end;

SetLength (V, n);

l:=0;

k:=0;

w:=0;

While w<=n-1 do begin

For i:=0 to m-1 do begin

r:=MNi[i];

For j:=0 to r-1 do begin

V[w]: =MXi[i];

w:=w+1;

end;

end;

end;

Sort;

end;

procedure OcenkaParam;

var i: integer;

alp1, alp2: real;

begin

X:=0;

SKv:=0;

for i:=0 to k-1 do begin

X:=X+(Xk[i]*OtnCh[i]);

end;

Form1. RzLabel27. Caption:=FloatToStrF (X, ffFixed, 10,4);

for i:=0 to k-1 do begin

SKv:=SKv+(sqr (Xk[i] - X));

end;

SKv:=SKv/(n-1);

S:=sqrt (SKv);

Form1. RzLabel25. Caption:=FloatToStrF (SKv, ffFixed, 10,4);

// Находим параметр альфа

alp1:=1+sqrt (1+(sqr (X)/SKv));

alp2:=1-sqrt (1+(sqr (X)/SKv));

If (alp1>0) then alp:=alp1

else begin

If (alp2>0) then alp:=alp2

else ShowMessage ('Параметр альфа — отрицательный');

end;

x0:=X*(alp-1)/alp;

Form1. RzLabel16. Caption:=FloatToStrF (x0, ffFixed, 10,4);

Form1. RzLabel17. Caption:=FloatToStrF (alp, ffFixed, 10,4);

If (alp>1) then begin

Mx:=alp*x0/(alp-1);

Form1. RzLabel18. Caption:=FloatToStrF (Mx, ffFixed, 10,4);

end

else Form1. RzLabel18. Caption:='Мат.ожидания нет';

If (alp>2) then begin

Dx:=alp*sqr (x0)/sqr (alp-1)*(alp-2);

Form1. RzLabel19. Caption:=FloatToStrF (Dx, ffFixed, 10,4);

end

else Form1. RzLabel19. Caption:='Дисперсии нет';

Form1. RzLabel15. Caption:=IntToStr (n);

end;

function F (x:real):real;

begin

F:=(1 — (Power ((x0/x), alp)));

end;

procedure Pareto;

var i: integer;

begin

If ((alp>0) and (x0>0)) then begin

For i:=0 to k-1 do begin

If (Xk[i]>x0) then begin

fx[i]: =(alp/x0) * (Power ((x0/Xk[i]), alp+1));

Frx[i]: =F (Xk[i]);

end

else begin

fx[i]: =0;

Frx[i]:=0;

end;

end;

end

else Form1. RzMemo2. Lines. Text:='Параметры распределения отрицательные';

// Делаем нормировку

Rx [k-1]: =Infinity;

Lx[0]:=x0;

If ((alp>0) and (x0>0)) then begin

For i:=0 to k-1 do begin

If (Lx[i]>x0) then begin

FLx[i]: =F (Lx[i]);

end

else FLx[i]: =0;

If (Rx[i]>x0) then begin

FRgx[i]: =F (Rx[i]);

end

else FRgx[i]: =0;

Px[i]:=FRgx[i] - FLx[i];

end;

end;

Sp:=0;

For i:=0 to k-1 do begin

Sp:=Sp+Px[i];

EmpCh[i]:=n*Px[i];

end;

Form1. RzLabel49. Caption:=FloatToStrF (Sp, ffFixed, 10,3);

For i:=1 to k do begin

Form1. RzStringGrid2. Cells [7, i]: =FloatToStrF (Px[i-1], ffFixed, 10,4);

Form1. RzStringGrid2. Cells [8, i]: =FloatToStrF (fx[i-1], ffFixed, 10,4);

Form1. RzStringGrid2. Cells [9, i]: =FloatToStrF (Frx[i-1], ffFixed, 10,4);

Form1. RzStringGrid2. Cells [11, i]: =FloatToStrF (EmpCh[i-1], ffFixed, 10,4);

end;

end;

procedure Gipotiza;

var i: integer;

begin

Form1. RzMemo2. Lines. Text:='';

Pareto;

HiNabl:=0;

For i:=0 to k-1 do begin

If (Px[i]<>0) then

HiNabl:=HiNabl+(Power ((Nk[i] - n*Px[i]), 2))/(n*Px[i]);

end;

Form1. RzLabel22. Caption:=FloatToStrF (HiNabl, ffFixed, 10,3);

Form1. RzLabel23. Caption:=FloatToStrF (HiKrit, ffFixed, 10,3);

a:= 'Гипотиза отброшена. Данная выборка не подчиняется распределению Парето';

b:= 'Гипотиза принята. Данная выборка подчиняется распределению Парето';

If (HiNabl>=HiKrit) then begin

Form1. RzMemo2. Text:=a;

Form1. RzLabel40. Caption:='>';

end

else begin

Form1. RzMemo2. Text:=b;

Form1. RzLabel40. Caption:='<';

end;

end;

procedure NewIntervals;

var i, m, l: integer;

t:real;

begin

If Form1. RzPageControl1. ActivePage=Form1. TabSheet1 then Vkladka1

else begin

If Form1. RzPageControl1. ActivePage=Form1. TabSheet2 then Vkladka2

else begin

If Form1. RzPageControl1. ActivePage=Form1. TabSheet3 then Vkladka3

end;

end;

Xmin:=V[0];

XMax:=V [n-1];

h:=(Xmax-Xmin)/k;

SetLength (Rx, k);

SetLength (Lx, k);

SetLength (Xk, k);

SetLength (Nk, k);

SetLength (OtnCh, k);

SetLength (Hk, k);

SetLength (Px, k);

SetLength (fx, k);

SetLength (Frx, k);

SetLength (FLx, k);

SetLength (FRgx, k);

SetLength (Fremp, k);

SetLength (EmpCh, k);

for i:=0 to k-1 do begin

Lx[i]: =Xmin+i*h;

Rx[i]:=Xmin+(i+1)*h;

Xk[i]:=(Lx[i]+Rx[i])/2;

end;

// подсчитываем число элементов, вошедших в интервал

m:=0;

l:=0;

Summa:=0;

While (m<=k-1) do begin

For i:=l to n-1 do begin

If (V[i]<=Rx[m]) then begin

Summa:=Summa+1;

l:=l+1;

end

else begin

Nk[m]: =Summa;

end;

Nk[m]:=Summa;

end;

Summa:=0;

m:=m+1;

end;

for i:=0 to k-1 do begin

OtnCh[i]: =Nk[i]/n;

Hk[i]:=OtnCh[i]/h;

end;

Sach:=0;

Sotch:=0;

Shk:=0;

for i:=0 to k-1 do begin

Sotch:=Sotch+OtnCh[i];

Shk:=Shk+h*Hk[i];

Sach:=Sach+Nk[i];

end;

Form1. RzLabel46. Caption:=FloatToStrF (Sach, ffFixed, 10,3);

Form1. RzLabel47. Caption:=FloatToStrF (Sotch, ffFixed, 10,3);

Form1. RzLabel48. Caption:=FloatToStrF (Shk, ffFixed, 10,3);

t:=0;

for i:=0 to k-1 do begin

Fremp[i]: =t+OtnCh[i];

t:=Fremp[i];

end;

For i:=1 to k do begin

Form1. RzStringGrid2. Cells [0, i]: =IntToStr (i);

Form1. RzStringGrid2. Cells [1, i]: =FloatToStrF (Lx[i-1], ffFixed, 10,3);

Form1. RzStringGrid2. Cells [2, i]: =FloatToStrF (Rx[i-1], ffFixed, 10,3);

Form1. RzStringGrid2. Cells [3, i]: =FloatToStrF (Xk[i-1], ffFixed, 10,3);

Form1. RzStringGrid2. Cells [4, i]: =IntToStr (Nk[i-1]);

Form1. RzStringGrid2. Cells [5, i]: =FloatToStrF (OtnCh[i-1], ffFixed, 10,4);

Form1. RzStringGrid2. Cells [6, i]: =FloatToStrF (Hk[i-1], ffFixed, 10,4);

Form1. RzStringGrid2. Cells [10, i]: =FloatToStrF (Fremp[i-1], ffFixed, 10,4);

end;

for i:=0 to k-1 do begin

Form1.diagr1. AddXY (Xk[i], Hk[i], FloatToStrF (Xk[i], ffFixed, 10,4), clBlue);

end;

OcenkaParam;

Gipotiza;

for i:=0 to k-1 do begin

Form1.diagr2. AddXY (Xk[i], fx[i], FloatToStrF (fx[i], ffFixed, 10,3), clRed);

end;