Использование структурных переменных

В простых вычислениях корректная подстановка аргументов не составляет проблемы, но при операциях с несколькими СВ и вложенными вызовами функций могут возникнуть технические трудности с передачей параметров распределений. Файл-функция f_Norm принимает параметры распределения m и  в полях m и s структурной переменной X:

function f=f_Norm(X,x)

f=f_Gauss((x-X.m)/X.s)/X.s;

Определив структурную переменную, ее можно передавать всем вызываемым функциям. Например, предыдущую команду можно изменить так:

>> X.m=2;X.s=3;x=-10:0.1:10; plot(x,f_Norm(X,x))

Класс нормально распределенных случайных величин

В классе объектов Norm_1 кроме параметров распределения можно определить и часто употребляемые функции обработки нормально распределенных СВ. В папке класса @Norm_1 содержатся следующие функции:

Norm_1(varargin) –	конструктор – создает объект класса с параметрами m, s, задаваемыми переменным списком (по умолчанию m = 0, s = 1);
setval(X,a,b) –	изменяет параметры объекта m, s;
Net(X,a,b,n,ns) –	разбивает интервал [a,b] на n равных частей (по умолчанию n=50, a, b вычисляются функцией Total, ns=4);
Total(X,a,b,ns) –	определяет границы интервала m  ns σ (в пределах [a,b], если они указаны);
f(X,x) –	вычислиет плотность распределения на сетке x;
Ver(X,a,b) –	вычисляет вероятность попадания СВ в интервал [a,b];
Gen(X,N) –	генерирует N случайных точек согласно закону распределения;
Fint(X,s,varargin) –	вычисляет полную вероятность с условным законом, заданным выражением в строке s, и параметрами в списке varargin;
display(X) –	используется средой MATLAB для вывода параметров объкта на экран

Создадим объект, соответствующий СВ X  N(3,2) и вычислим вероятность отклонения от МО не более, чем E в большую сторону, по определению (равную 0,25):

>> X=Norm_1(3,2); p=Ver(X,3,3+2*sqrt(2)*0.477)

p = 0.2500

Генерируем 10000 случайных точек и вычислим частоту попаданий в тот же интервал:

>> N=10000;t=Gen(X,N);p=sum(t>3 & t<3+2*sqrt(2)*0.477)/N

p = 0.2502

Литература к лекции

1. Ришняк А.Г. Овчинников А.Ф. Вероятностные задачи теории эффективности действия: Учеб. пособие. – М: Изд-во МГТУ им. Баумана, 2006. 304 с.

Контрольные вопросы и задачи

1. Какому закону подчиняется расстояние между ближайшими точками простейшего пуассоновского поля на плоскости?

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

3. Какова характерная особенность числовых характеристик показательного закона?

4. Как определяется надежность последовательного и параллельного соединений элементов по характеристикам их надежности? Какому закону подчиняются потоки отказов в таких соединениях?

5. В чем суть объектно-ориентированного моделирования надежности сложных соединений элементов?

6. Объясните содержание электроннй формулы MassDyn для моделирования занятости многоканальной системы массового обслуживания. Как можно анализировать вероятность немедленной обработки заявки в системах массового обслуживания с очередями и отказами?

7. Назовите причины, по которым неизвестное распределение может быть принято нормальным.

ПРИЛОЖЕНИЕ к лекции 5

Листинг 5.1. Функция P_Gauss для вычисления функции распределения стандартного нормального закона:

%Вычисляет функцию Гаусса

function F = P_Gauss(x)

F=F_LaplasV(x)+0.5;

Листинг 5.2. Функция f_Norm(X,x) вычисляет плотность нормального распределения:

% Вычисляет плотность нормального распределения с параметрами X.m, /X.s

% в точках, заданных вектором x.

function f=f_Norm(X,x)

f=f_Gauss((x-X.m)/X.s)/X.s;

Листинг 5.3. Функция MassObs (X,x) вычисляет время ожидания в очереди в многоканальной системе массового обслуживания:

function [D,A]=MassObs(L,T,N,n)

if nargin<4 n=1; end

A=random('exp',L,N,1); A=cumsum(A);

B=random('exp',T,N,1);

D=zeros(size(A));T=zeros(1,n);

for i=1:N

[Tr, I]=min(T);

if A(i) > Tr

T(I)=A(i)+B(i);

else

D(i)=Tr-A(i);

T(I)=Tr+B(i);

end

end

Листинг 5.4. Функции класса надежности Rab.

Конструктор класса Rab:

function X=Rab(varargin)

X=class(struct('Class','Rab','P',[],'S',''),'Rab');

X=setval(X,varargin{:});

Конструктор класса Rab:

function X=setval(X,varargin)

for i=1:nargin-1

x=varargin{i};

if ischar(x) X.S=x;

elseif isa(x,'Rab')

X=Ap(X,x,50000);

else

if length(x)==1 x=[x 1 1]; end

X.P(end+1,:)=x;

end

end

function X=Ap(X,A,N)

[x,f,F]=Approx(A,N);

Goal=inline('sum((p_Weibull(x,L(1),L(2))-F).^2)','L','x','F');

d=diff(x([2 3]));

L=fminsearch(Goal,[1,1],optimset('MaxIter',50,'TolX',0.0001), [x x(end)+d],F);

X=Rab([L 1]);

Конструктор класса Rab:

function out=CalcP(X,P)

V=RabLog(P);

k=1;

for i=1:length(X.P)

m=X.P(i,3);

U(i)=sum(V(k:k+m-1));

k=k+m;

end

out=Value(eval(MkStr(X.S)));

%

function S=MkStr(s)

a=(s>='0' & s<='9');

b=diff([0 a]);

beg=find(b==1);

kon=find(b==-1);

m=length(beg); beg(m+1)=length(s)+1;

if length(kon)<m

kon(m)=length(b)+1;

end

S=s(1:beg(1)-1);

for i=1:m

S=strcat(S, ['U(' s(beg(i):kon(i)-1) ')' s(kon(i):beg(i+1)-1)]);

end

Функция Gen генерирует N реализаций мометов отказов элементов системы:

function [V,M]=Gen(X,N)

if nargin<2 N=1;end

n=length(X.P);

M = zeros(N,n);V=zeros(1,N);k=1;

for i=1:n

m=X.P(i,3);

C=Gen('Wei',X.P(i,1:2),N,m);

M(:,k:k+m-1)=C;k=k+m;

end

for i=1:N

V(i)=CalcP(X,M(i,:));

end

Функция Approx строит гистограмму и эмпирическую функцию распределения по создаваемым ею N реализациям мометов отказов элементов системы:

function [x,f,F]=Approx(X,N)

if nargin<2 N=10000; end

P=sort(Gen(X,N));

u=P(end);

d=u/30;x=0:d:u;f=histc(P,x)/N; F=[0 cumsum(f)];f=f/d;

Функция MassModel моделирует вычисляет продолжительность ожидания в очереди (queue = 1) или вероятность обработки заявки (queue = 0) в n-канальной СМО:

function [D,A]=MassModel(law1,law2,N,n,queue)

if nargin<5 queue=1; end

if nargin<4 n=1; end

A=Gen(law1{:},N,1); A=cumsum(A);

B=Gen(law2{:},N,1);

D=zeros(size(A));T=zeros(1,n);

for i=1:N

[Tr, I]=min(T);

if A(i) > Tr

T(I)=A(i)+B(i);

else

if queue==0

D(i)=1;

else

D(i)=Tr-A(i);

T(I)=Tr+B(i);

end

end

end

Функция MassDyn моделирует занятость n-канальной СМО с отказами (queue = 0) и очередями (queue = 1):

function Q=MassDyn(law1,law2,N,n,queue,T)

Q=zeros(1,length(T)-1);K=Q;

for k=1:N

[a,b]=MassDyn_(law1,law2,500,n,queue,T);

Q=Q+a;

K=K+b;

end

if queue==1 Q=Q./K;else Q=(1-Q)./K;end

function [Q,K,D,A]=MassDyn_(law1,law2,N,n,queue,T)

[D,A]=MassModel(law1,law2,N,n,queue);

nt=length(T);

Q=zeros(1,nt-1);K=Q;

for i=1:nt-1

I=find(A>T(i) & A<T(i+1));

K(i)=length(I);

Q(i)=sum(D(I));

end

Листинг 5.5. Функции класса Norm_1.

Конструктор класса Norm_1.

function X=Norm_1(varargin)

X=class(struct('Class','Norm_1','m',0,'s',1),'Norm_1');

X=setval(X,varargin{:}); end

Функция Norm_1/setval изменяет параметры объекта:

function X=setval(X,a,b)

if nargin>1 & isstruct(a)

X.m=a.m; X.s=a.s;

else

if nargin<3 b=1; end

if nargin<2 m=0; end

X.m=a; X.s=b;

end

Функция Norm_1/ f уточняет границы интервала вероятных значений:

function out=f(X,x)

out=f_Gauss((x-X.m)/X.s)/X.s;

Функция Norm_1/Ver вычисляет вероятность попадания в интервал:

function P=Ver(X,a,b)

P=p_Norm(struct(X),a,b);

%

function P=p_Norm(X,x1,x2)

m=[prod(size(x1)),prod(size(x2)),prod(size(X.m)),prod(size(X.s))];

[n, I]=max(m);

switch I

case 1, S=zize(x1); break

case 2, S=size(x2);

case 3, S=size(X.m);

case 4, S=size(X.s);

end

P=zeros(S);

for i = 1:n

if i<= m(1) x1i=x1(i); end

if i<= m(2) x2i=x2(i); end

t=linspace(x1i,x2i,100);

if i<= m(3) Y.m=X.m(i); end

if i<= m(4) Y.s=X.s(i); end

p=fNorm(Y,t);

P(i)=defint(t,p);

end

Функция Norm_1/ Total уточняет границы интервала вероятных значений:

function [a, b]=Total(X,a,b)

a0=X.m-4*X.s; b0=X.m+4*X.s;

if (nargin<2 | a<a0) a=a0; end

if (nargin<3 | b>b0) b=b0; end

Функция Norm_1/Fint вычисляет полную вероятность с условным законом, заданным выражением в строке s, и параметрами в списке varargin:

function out=Fint(X,str,varargin)

[a, b]=full(X,varargin{:});

x=linspace(a,b,50);

out=Trap(f(X,x).*eval(str),x);

Листинг 5.8. Функция Norm_1/display используется средой MATLAB для вывода параметров объкта на экран:

function display(X)

line1=[' Norm_1' sprintf(' %s', inputname(1))];

line2= sprintf('MO = %0.3f',X.m);

line3= sprintf(' CKO = %0.3f',X.s);

disp([line1 ' [' line2 line3 ']' ])

Листинг 5.3. Функция Norm_1/Gen генерирует N случайных точек:

function X=Gen(X,N)

if nargin<2 N=1; end

X=randn(1,N).*X.s + X.m;