Объектная среда для вероятностного анализа в системе MATLAB базируется на понятии "электронная формула". Электронные формулы в рабочей области MATLAB синтаксически оформляются как вызовы функций, имеющих содержательное толкование в контексте теории вероятностей и выполняющие действия, зависящие от комбинации аргументов. Аргументы электронных формул, как правило, различаются по смыслу, будучи объектами определенных классов: событий, случайных величин, геометрических фигур и т.д. Вычисления при исполнении функций выполняют объекты во взаимодействии между собой. Математический аппарат теории вероятностей реализован иерархией классов, исходные данные вероятностных задач принимаются конструкторами классов, создающих объекты в контексте задачи, электронные формулы управляют взаимодействием объектов и формируют результат, соответствующий комбинации аргументов.

Включает в себя следующие разделы:

• вероятности сложных событий;
• формула Байеса;
• формула Бернулли;
• дискретные случайные величины;
• базовый класс непрерывных случайных величин;
• показательный закон и связанные с ним классы распределений;
• нормальный закон распределения;
• класс двумерных нормально распределенных векторов;
• функции случайных величин;
• мультипликативные функции случайных величин;
• аддитивные функции случайных величин.

MATLAB версии не ниже R2019a

Выполните в MATLAB

eval(webread('https://exponenta.ru/install/m4'))

Скачайте и запустите файл ProbabilityAnalysisM4.mltbx

Перед использованием функций библиотеки выполните команду addVerToPathShort

Выполните в MATLAB

m4update.m

Если X – случайная величина с известным законом распределения, G – геометрическая фигура или массив непересекающихся фигур, электронная формула [P, p] = Ver(X, G) вычисляет вероятность попадания P случайной величины X в область G. Объект X создается конструктором одного из классов случайных величин с соответствующими параметрами, например, X = Norm_2([1;2],[3,4],0.5) – это объект нормального распределения на плоскости с математическими ожиданиями по координатам, равными 1 и 2, среднеквадратическими отклонениями 3 и 4, коэффициентом корреляции между координатами 0,5. Объект G может быть создан конструктором кругов заданного радиуса G = Circ(1), прямоугольников G = Rect(2,3) или других фигур, может быть также результатом объединения или пересечения фигур. Результат выполнения электронной формулы Ver зависит от структуры объекта G. Если G – массив фигур, P – вероятность попадания хотя бы в одну из них, p – массив вероятностей попадания в отдельные фигуры. Из массива вероятностей p конструктор класса случайных событий Randevent создает массив событий с вероятностями из массива p: A = Randevent(p). Выражение для суммы событий sum(A) в данном примере даст результат, совпадающий со значением P, так как элементы массива A – несовместные события.

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

>> [A,B,C]= Randevent(0.3, 0.4, 0.5); D=A+B*C, E=(A+B)*C
Событие D: P(D) =   0.44
Событие E: P(E) =   0.29

В операциях умножения учитывается зависимость между событиями, заданная условными вероятностями P(A|C) = 0.45, P(B|C) = 0.35:

>> A=Set(A,C,0.45); B=Set(B,C,0.35); D=A+B*C, E=(A+B)*C
Событие D: P(D) =  0.396
Событие E: P(E) =  0.321

В операциях с массивами событий функция sum (сложение) учитывает совместность слагаемых, функция prod (произведение) – попарную зависимость:

>> p=[0.2 0.4 0.45 0.25 0.5 0.1 0.6 0.15 0.4]; 
>> SYS = Randevent(p), A = sum(SYS),B=prod(SYS(1:4))+prod(SYS(5:9))
Группа событий SYS: P(SYS(1:9)) =    0.2   0.4  0.45  0.25   0.5   0.1   0.6  0.15   0.4
Событие A: P(A) =  0.982
Событие B: P(B) = 0.0108

Если операндами в формуле события U служат индексы элементов массива событий '1 + 2 + (3+4)*(5+6) + (7+8)*9', объект класса Randevent формирует алгебраическое выражение из элементов массива событий согласно формуле из индексов, при вычислении которого используются переопределенные двухместные операции:

>> U = Set(SYS, '1+2+(3+4)*(5+6)+(7+8)*9')
Событие U: P(U) =  0.761

Электронная формула Bayes класса Randevent вычисляет апостериорную вероятность интересующей гипотезы по последовательности событий (A, A, ..., Ā, ...), наблюдавшихся в контрольных испытаниях. В качестве первого аргумента она получает объект уточняемой гипотезы Hi и последовательность объектов A и not(A), возвращаемый результат – апостериорная вероятность уточняемой гипотезы:

HiPost=Bayes(Hi, A, ..., not(A), ...).

Порядок чередования событий и их отрицаний в списке аргументов безразличен, можно заменить группу одинаковых событий парой аргументов – (событие, кратность):

HiPost=Bayes(Hi, A, n, not(A), m).

Вычисления по частной, общей и обобщенной формулам Бернулли выполняет электронная формула Ber: out = Ber(p, n, m), где p – вероятность успеха в каждом испытании, или вектор-строка вероятностей успеха или вектор-столбец вероятностей каждого из множества возможных исходов; n – число испытаний; m – интересующее число успехов, или вектор-строка диапазона интересующего числа успехов или вектор-столбец интересующей комбинации чисел исходов. Выходная величина out в прямых задачах имеет значение вероятности наступления m событий, в обратных задачах – принимает значение вычисленной величины, заданной пустым вектором [ ] в списке аргументов. Электронная формула Ber с пустым первым аргументом и заданной четвертым аргументом доверительной вероятностью Dov вычисляет доверительный интервал [p1, p2] для вероятности успеха по m успешным исходам в n испытаниях.

В объектно-ориентированной вероятностной модели объекты распределения участвуют в алгебраических операциях, вычисляют числовые характеристики и вероятности связанных с ними событий. Теоретические законы распределения, различающиеся набором параметров и их интерпретацией, реализованы в производных классах от DRV. Конструктор DRV без аргументов выводит сведения об основных функциях класса и список производных классов:

>> DRV
 конструкторы дискретного распределения, создающие объект X:
  X=DRV(x;p) или X=DRV([x;p])-  x - массив значений; p - массив вероятностей
  X=DRV(A); A-случайное событие (объект класса Randevent), X-индикатор события A
 основные функции класса:
  Val(X,Ind) - возможные значения СВ X (с индексами Ind, если 2-й аргумент задан)
  Ver(X,Ind) - вероятности значений СВ X (с индексами Ind, если 2-й аргумент задан)
  MO(X), SKO(X) - вывод МО и СКО СВ X)
  Show(X1,s1,X2,s2,...) - вывод многоугольников распределения СВ X1,X2,...)
 конструкторы производных классов:
  BIN(p,n) - биномиальное распределение
  POI( L ) - распределение Пуассона
  GEO(p,d) - геометрическое распределение (d=0) или сдвинутое геометрическое (d>0)
  HyperGEO(N,R,L) - гипергеометрическое (R дефектных из N, L - объем выборки)

Массив вероятностей, принимаемых конструктором, проверяется на нормированность единицей. При этом, если один из элементов массива вероятностей задан числом -1, он заменяется дополнением суммы остальных элементов до единицы. Если модуль отрицательного числа меньше единицы, элемент заменяется абсолютным значением, а дополнение до единицы равномерно распределяется между нулевыми элементами массива. Аргументом конструктора DRV может быть случайное событие, тогда объект создается как индикатор этого события. Алгебраические операции с объектами СВ выполняются по правилам теории вероятностей, сохраняют свойства распределений, такие как устойчивость закона Пуассона к сложению:

>> X2=POI(2); X3=POI(3); X5=X2+X3
X5: Распределение Пуассона (5) MO=5, SKO=2.24

Базовый класс CRV непрерывных случайных величин содержит закон распределения в виде функции плотности вероятности или функции распределения, или их дискретных значений на сетке. Базовый класс порождает классы основных законов распределения:

>> CRV
 конструкторы непрерывного распределения, создающие объект X:
  X=CRV(x;p) или X=DRV([x;p])-  x - массив значений; p - массив вероятностей
  X=DRV(A); A-случайное событие (объект класса Randevent), X-индикатор события A
 основные функции класса:
  Val(X,Ind) - возможные значения СВ X (с индексами Ind, если 2-й аргумент задан)
  Ver(X,Ind) - вероятности значений СВ X (с индексами Ind, если 2-й аргумент задан)
  MO(X), SKO(X) - вывод МО и СКО СВ X)
  Show(X1,s1,X2,s2,...) - вывод многоугольников распределения СВ X1,X2,...)
  EXP( L )   - показательное распределение
  ERL( L,k ) - распределение Эрланга порядка k
  WEI(p1,p2) - распределение Вейбулла
  RAYL(p)    - распределение Рэлея
  Norm_1(m,s)- нормальное распределение
  RND(a,b)   - равномерное распределение
  GAMMA(a,b) - Гамма-распределение

Базовый класс CRV создает объекты произвольных распределений, заданных формулой или статистическим рядом. Например, случайная длина проекции отрезка длиной l может быть представлена объектом класса CRV, заданным интегральной функцией F(x) = 2/pi*asin(x/L) или плотностью распределения f(x) = 2/pi/sqrt(L^2 – x^2):

>> L = 10; XF=CRV('2/pi*asin(x/L):L',[0,L]), Xf=CRV('2/pi./sqrt(L^2-x.^2):L',[0,L])
XF: Распределение F:2/pi*asin(x/10) [0 10]: MO=6.36, SKO=3.03
Xf: Распределение f:2/pi./sqrt(100-x.^2) [0 10]: MO=6.36, SKO=3.08

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

Вероятностный анализ событий, связанных с пуассоновскими потоками отказов, заявок, сигналов может быть выполнен без методических упрощений объектами класса показательного закона EXP и родственных классов закона Эрланга ERL и Вейбулла WEI. В следующем примере конструктор EXP создает объект показательного распределения с параметром 2, объекты выполняют алгебраические действия по правилам теории вероятностей, функция Show строит графики плотности и функции распределения:

>> T=EXP(1.5), Show(T,'Fk', T, 'fr'), T3=T+T+T, figure, Show(T,'k',T+T,'r',T3,'b',T3+T,'g')
T: Распределение экспоненциальное (1.5)  [0 Inf]: MO=0.67, SKO=0.67
T3: Распределение Эрланга порядка 2(1.5) MO=2, SKO=1.15

В классах непрерывных СВ перемножение функций распределения выполняет программа Prodf, создающая объект с сеточным распределением класса CRV. Следующая команда определяет массив объектов X класса EXP массивом параметров, аппроксимирует произведение Prodf(X) = Y объектом класса WEI, выводит графики функций распределения объектов X, плотности и функции распределения СВ Y, а также функцию распределения Вейбулла (точками), близкую к графику функции распределения СВ Y:

>> X=EXP(1.2:0.2:1.8), Y=Prodf(X), W=WEI(Y), Show(X,'Fb', Y, 'r',Y,'Fk',W,'Fr.')
Массив X:
 Распределение экспоненциальное (1.2)  [0 Inf]: MO=0.83, SKO=0.83
 Распределение экспоненциальное (1.4)  [0 Inf]: MO=0.71, SKO=0.71
 Распределение экспоненциальное (1.6)  [0 Inf]: MO=0.63, SKO=0.63
 Распределение экспоненциальное (1.8)  [0 Inf]: MO=0.56, SKO=0.56
Y: Распределение  [0, Inf]: МО =1.37, SKO = 0.83
W: Распределение Вейбулла (0.6416,1.7807) MO=1.39, SKO=0.81

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

Универсальную модель занятости многоканальной системы массового обслуживания (СМО) с отказами и очередями воспроизводит стохастическая процедура, реализованная в программе Qsystem. Она получает объекты потоков заявок X, обслуживания Y и число каналов n. Признаком СМО с ожиданием служит отрицательный знак аргумента n. Дополнительными необязательными аргументами задаются число повторений процесса N для набора статистики (по умолчанию, 10000), длительность процесса T (10) и число net (50) моментов времени, в которые оценивается вероятность наличия свободного канала. Программа N раз имитирует последовательность поступления заявок, длительность их обслуживания, отказы в обслуживании или постановку в очередь. По накопленной статистике восстанавливается зависимость незанятого состояния СМО с отказами или среднего времени ожидания в СМО с очередями. Для аппроксимации зависимости подбирается подходящая функция, которая выводится в качестве результата вместе с параметрами. В одноканальной СМО с очередью при плотностях входного и выходного потоков λ = 1, μ = 1/0,99 = 1.01 время ожидания стабилизируется на уровне 98, но установившийся режим наступает примерно через 100000 единиц времени, при не менее 10000 повторений испытаний (в примере - 1000 испытаний для ускорения расчета):

>> X=EXP(1); Y=EXP(1.01);T=logspace(-1,5,50);[Fn,Par]=Qsystem(X,Y,-1,1000,T);t=Par(2),plot(T,Fn(T,Par))
t  =  70.1634

Сведения об основных функциях класса NORM можно получить по имени класса.

>> NORM
 X=NORM(m,s) - конструктор объектов нормального распределения
 X=NORM(m) или X=NORM - конструктор объектов (по умолчанию m=0, s=1)
   основные функции класса:
 y = f(X,x) - вычисляет плотность y распределения СВ X на сетке x
 x = Net(X,n,N) - строит сетку n точек на интервале "N сигм" (по умолчанию n=50,N=4)
 T = Gen(X,m,n,Y) - генерирует nxm-матрицу T реализаций X, смещенных на Y (n=1,Y=[])
 F=fint(X,str) - вычисление по полной вероятности str - выражение условной вероятности

Конструктор Norm_2 создает объект, выполняющий все действия, предусмотренные математической моделью двумерных нормальных распределений. Класс Norm_2 наследует базовый класс геометрических объектов вместе с методами пространственных преобразований (сдвиг, масштабирование, повороты). Вектор математического ожидания (МО) как центр рассеивания позиционирует объект, среднеквадратическое отклонение (СКО) определяет размеры единичного эллипса рассеивания. Двумерный вектор-столбец принимается конструктором как вектор МО, положительный [1, 2]-вектор – как СКО. Скалярный аргумент учитывается конструктором как коэффициент корреляции. Сведения о способах определения объектов и методах работы с ними можно получить вызовом конструктора Norm_2 без аргументов и принимающей переменной:

>> Norm_2
   Конструкторы объектов нормального распределения:
 X=Norm_2(m,K); m - 2x1-вектор МО =[0;0], K – корреляционная матрица
 X=Norm_2(m,s,r),  m - 2x1-вектор МО =[0;0], s - 1x2-вектор СКО =[1;1], r - коэф.корреляции =0
 X=Norm_2(X1,X2,r), X1, X2 - объекты класса NORM, r - коэф.корреляции =0
   Основные функции класса:
 dens = f(X, x, y) - вычисляет плотность dens распределения СВ X на сетке x, y
 [x,y]=Net(X,n,N) - строит сетку n точек на интервале "N сигм" (по умолчанию n=50,N=4)
 T = Gen(X,m,n) - генерирует nxm-матрицу T реализаций X
 Y=Rot(X,a1,a2) - поворот осей на угол a1[град] (или a2[рад], если a1=[])
 P=Ver( X, G) – вычисляет вероятность P попадания случайной точки X в область G

Переменную в левой части конструктор без параметров определяет объектом стандартного распределения с единичными СКО без систематических ошибок и корреляции. Операции сдвига, масштабирования и поворота иллюстрируют операции с объектом Y0:

>> Y0=Norm_2, Y1 = Y0*[2, 1], Y2 = Y1+[12;3], Y3 = Rot(Y1,-30) + [12; 3]
 Norm_2 Y0: MO = [0 0], CKO = [1 1]
 Norm_2 Y1: MO = [0 0], CKO = [2 1]
 Norm_2 Y2: MO = [12 3], CKO = [2 1]
 Norm_2 Y3: MO = [12 3], CKO = [1.8028 1.3229], r = 0.54

Электронная формула Ver(X, G) вычисляет вероятность попадания нормально распределенной точки X в геометрическую область G численным интегрированием или по специальным формулам в частных случаях.

Аналитические преобразования функции случайных величин могут быть заменены универсальной процедурой, обрабатывающей интерпретируемое выражение функции согласно теоретической формуле. Электронная формула FUN получает выражение функции S, один или несколько случайных аргументов X1, ..., параметры a1, ..., и создает объект распределения Y класса CRV: Y=FUN(S, X1, ..., a1, ...) Выражение S для площади проекции параллелепипеда с параметрами a, b, c содержит перечень параметров в том порядке, в каком они принимаются из списка входных данных, угол Fi, равномерно распределенный определен объектом класса равномерного распределения RND, объект распределения угла Teta – по закону косинуса:

>> S = 'b*c*sin(Teta) + a*(c*abs(sin(Fi)) + b*abs(cos(Fi))).*cos(Teta) : a,b,c'; a=9; b=6; c=4;
>> Teta = CRV('cos(x)', [0, pi/2]); Fi = RND([0, 2*pi]); Pr3 = FUN(S, Teta, Fi, a, b, c)
Pr3: Распределение  [24.003, 69.1953]: МО =57.05, SKO = 9.39

Для вычисления площади проекции S геометрического объекта G класса многогранников PolyGran на плоскость, ориентированную под сферическими углами Teta, Fi, можно использовать электронную формулу ProectS: S = ProectS(G, Teta, Fi)

Мультипликативные операции с объектами распределений учитывают особенности сомножителей. В примере построенные законы распределения площади прямоугольника, длины сторон которого X1 ∈ N(10, 2) и X2 ∈ N(15, 4) в одном случае независимы, а двух других зависимы с коэффициентами корреляции 0,5 и –0.5:

>> r=     0; X=Norm_2([10;15],[2 4], r); [X1,X2]=X12(X); Y1=X1*X2; Show(Y1, 'b')
>> r=  0.5; X=Norm_2([10;15],[2 4], r); [X1,X2]=X12(X); Y2=X1*X2; Show(Y2, 'r')
>> r= -0.5; X=Norm_2([10;15],[2 4], r); [X1,X2]=X12(X); Y3=X1*X2; Show(Y3, 'k')

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

>> R1=RND(0,1),R2=R1+R1, R3=R2+R1,R4=R3+R1, R6=R4+R2, N=NORM(R6), Show(R1,R2,R3,R4,R6,N ,'.')
R1: Распределение RND [0 1]: MO=0.5, SKO=0.29
R2: Распределение Симпсона[0, 2]: МО =1, SKO = 0.41
R3: Распределение  [0, 3]: МО =1.5, SKO = 0.5
R4: Распределение  [0, 4]: МО =2, SKO = 0.58
R6: Распределение  [0, 6]: МО =3, SKO = 0.71
N: Распределение нормальное (3 0.71) : MO=3, SKO=0.71

Композиция нормального и равномерного законов выполняется операцией сложения с объектами классов NORM и RND:

P = Ver(NORM(m, sigma) + RND(a,b), [z1, z2])