7. Додаток

У стандартному Statistics Toolbox системи MATLAB немає програм для обчислення функції розподілу, густини розподілу і квантилів розподілу Колмогорова. Нижче приведені процедури для обчислення цих функцій. Їх потрібно скопіювати (наприклад, через Clipboard) у файли з відповідними іменами, і помістити ці файли в будь-який каталог, доступний системі MATLAB.

Щоб перевірити, чи доступна функція, наберіть в MATLAB Command Window команду

help [ім'я функції]

Якщо функція доступна, на екран буде видана довідкова інформація по функції.

Текст файлу kolmcdf.m (функція розподілу Колмогорова).

function p = kolmcdf(x)

% KOLMCDF Kolmogorov cumulative distribution function (cdf).

% P = KOLMCDF(X) computes the Kolmogorov cdf at the values in X.

% The size of P is the common size of X. A scalar input

% functions as a constant matrix of the same size as the other inputs.

if nargin < 1,

error('Requires at least one input argument.');

end

p=zeros(size(x));

num=find(x>0);

xnum=x(num);

pnum=ones(size(num));

for k=1:1000,

add=(-1)^k*exp(-2*k^2*xnum.^2);

pnum=pnum+2*add;

if norm(add,1)==0,

break

end

end

p(num)=p(num)+pnum;

Текст файлу kolmpdf.m (густина розподілу Колмогорова).

function p = kolmpdf(x)

% KOLMPDF Kolmogorov probability distribution function (pdf).

% P = KOLMPDF(X) computes the Kolmogorov pdf at the values in X.

% The size of P is the common size of X. A scalar input

% functions as a constant matrix of the same size as the other inputs.

if nargin < 1,

error('Requires at least one input argument.');

end

p=zeros(size(x));

num=find(x>0);

xnum=x(num);

pnum=zeros(size(num));

for k=1:1000,

add=4*(-1)^(k+1)*k^2*exp(-2*k^2*xnum.^2).*xnum;

pnum=pnum+2*add;

if norm(add,1)==0,

break

end

end

p(num)=p(num)+pnum;

Текст файлу kolminv.m (квантилі розподілу Колмогорова).

function x = kolminv(p)

% KOLMINV Inverse of the Kolmogorov cumulative distribution function.

% P = KOLMINV(X) returns the inverse of the Kolmogorov

% cumulative distribution function, at the values in P.

% The size of X is the common size of P. A scalar input

% functions as a constant matrix of the same size as the other input.

if nargin < 1,

error('Requires at least one input argument.');

end

% Initialize X to zero.

x=zeros(size(p));

% Return NaN if the arguments are outside their respective limits.

k = find(p < 0 | p > 1);

if any(k),

tmp = NaN;

x(k) = tmp(ones(size(k)));

end

k = find(p == 0);

if any(k),

x(k) = zeros(size(k));

end

k = find(p==1);

if any(k),

tmp = inf;

x(k) = tmp(ones(size(k)));

end

% Newton's Method.

% Permit no more than count_limit interations.

count_limit = 50;

count = 0;

k = find(p > 0 & p < 1);

If isempty(k)

return;

end

pk = p(k);

% Use the mean as a starting guess.

xk = (pk-0.05)/0.9+0.5;

h = ones(size(pk));

crit = sqrt(eps);

% Break out of the iteration loop for the following:

% 1) The last update is very small (compared to x).

% 2) The last update is very small (compared to 100*eps).

% 3) There are more than 50 iterations. This should NEVER happen.

while(any(abs(h) > crit * abs(xk)) & max(abs(h)) > crit ...

& count < count_limit),

count = count+1;

h = (kolmcdf(xk) - pk) ./ kolmpdf(xk);

xnew = xk - h;

xk = xnew;

end

% Return the converged value(s).

x(k) = xk;

if count==count_limit,

fprintf('\nWarning: KOLMINV did not converge.\n');

str = 'The last step was: ';

outstr = sprintf([str,'%13.8f'],h);

fprintf(outstr);

end

12