In [43]:

    

%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
import sympy as sp
import sympy.printing.mathml as mathml
from sympy.printing.mathml import print_mathml

sp.init_printing(use_unicode=True)


    

In [20]:

    

print('Задача 1. Моделирование неоднородной дискретной цепи Маркова.')
k = sp.Symbol('k')
print('P(k)=')
P_k = sp.Matrix([
        [(k + 3) / (2 * (k + 2)), (k-1) / (2 * (k + 2)), 1 / (k + 2)],
        [1 / (k + 1), (k - 1) / (k + 1), 1 / (k + 1)],
        [2 / (k + 2), 1 / (k + 2), (k - 1) / (k + 1)],
])
fP_k = lambda a: P_k.replace(k, a)

P = sp.Matrix([[0.4, 0.2, 0.4], [0.3, 0.5, 0.2], [0.2, 0.3, 0.5]])

P_k, P


    

    




Задача 1. Моделирование неоднородной дискретной цепи Маркова.
P(k)=






    
Out[20]:





$$\left ( \left[\begin{matrix}\frac{k + 3}{2 k + 4} & \frac{k - 1}{2 k + 4} & \frac{1}{k + 2}\\\frac{1}{k + 1} & \frac{k - 1}{k + 1} & \frac{1}{k + 1}\\\frac{2}{k + 2} & \frac{1}{k + 2} & \frac{k - 1}{k + 1}\end{matrix}\right], \quad \left[\begin{matrix}0.4 & 0.2 & 0.4\\0.3 & 0.5 & 0.2\\0.2 & 0.3 & 0.5\end{matrix}\right]\right )$$

In [25]:

    

# fP_k = sp.lambdify(k, P_k)
P_k.multiply(P)


    

    
Out[25]:





$$\left[\begin{matrix}\frac{0.3 k - 0.3}{2 k + 4} + \frac{0.4 k + 1.2}{2 k + 4} + \frac{0.2}{k + 2} & \frac{0.5 k - 0.5}{2 k + 4} + \frac{0.2 k + 0.6}{2 k + 4} + \frac{0.3}{k + 2} & \frac{0.2 k - 0.2}{2 k + 4} + \frac{0.4 k + 1.2}{2 k + 4} + \frac{0.5}{k + 2}\\\frac{0.3 k - 0.3}{k + 1} + \frac{0.6}{k + 1} & \frac{0.5 k - 0.5}{k + 1} + \frac{0.5}{k + 1} & \frac{0.2 k - 0.2}{k + 1} + \frac{0.9}{k + 1}\\\frac{0.2 k - 0.2}{k + 1} + \frac{1.1}{k + 2} & \frac{0.3 k - 0.3}{k + 1} + \frac{0.9}{k + 2} & \frac{0.5 k - 0.5}{k + 1} + \frac{1.0}{k + 2}\end{matrix}\right]$$

In [47]:

    

sp.printing.latex(_25)
# print('<?xml version="1.0" encoding="UTF-8"?><math xmlns="http://www.w3.org/1998/Math/MathML">' + print_mathml(_25) + "</math>")


    

    
Out[47]:





'\\left[\\begin{matrix}\\frac{0.3 k - 0.3}{2 k + 4} + \\frac{0.4 k + 1.2}{2 k + 4} + \\frac{0.2}{k + 2} & \\frac{0.5 k - 0.5}{2 k + 4} + \\frac{0.2 k + 0.6}{2 k + 4} + \\frac{0.3}{k + 2} & \\frac{0.2 k - 0.2}{2 k + 4} + \\frac{0.4 k + 1.2}{2 k + 4} + \\frac{0.5}{k + 2}\\\\\\frac{0.3 k - 0.3}{k + 1} + \\frac{0.6}{k + 1} & \\frac{0.5 k - 0.5}{k + 1} + \\frac{0.5}{k + 1} & \\frac{0.2 k - 0.2}{k + 1} + \\frac{0.9}{k + 1}\\\\\\frac{0.2 k - 0.2}{k + 1} + \\frac{1.1}{k + 2} & \\frac{0.3 k - 0.3}{k + 1} + \\frac{0.9}{k + 2} & \\frac{0.5 k - 0.5}{k + 1} + \\frac{1.0}{k + 2}\\end{matrix}\\right]'

In [ ]:

    

for n = 1:10
    sost = 1;
for i = 1:20
  P = [(i+3)/(2*(i+2)) (i-1)/(2*(i+2)) 1/(i+2);
        1/(i+1) (i-1)/(i+1) 1/(i+1);
        2/(i+2) 1/(i+2) (i-1)/(i+1)];      
  
  alpha1(i) = rand(1); 
  
  if(alpha1(i) < P(sost,1))
    Res(n,i) = 1;
    sost = 1;
  else if((P(sost,1) < alpha1(i)) && (alpha1(i) <= P(sost,1) + P(sost,2)))
        Res(n,i) = 2;
        sost = 2;
  else if(alpha1(i) > P(sost,1) + P(sost,2))
        Res(n,i) = 3;
        sost = 3;
  end
  end
  end
end
end

disp('');
disp('Результат: ');
disp(Res);


disp('-----------------------------------------');
disp('Задача 2. Моделирование однородной дискретной цепи Маркова.');
disp(' ');
disp('P=');
P = [0.1 0.3 0.6;
        0.7 0.1 0.2;
        0.1 0.8 0.1];
disp(P);    
disp(' ');


for n = 1:10
    
    sost = 1;
for i = 1:20    
  
  alpha2(i) = rand(1); 
  
  if(alpha2(i) < P(sost,1))
    Res(n,i) = 1;
    sost = 1;
  else if((P(sost,1) < alpha2(i)) && (alpha2(i) <= P(sost,1) + P(sost,2)))
        Res(n,i) = 2;
        sost = 2;
  else if(alpha2(i) > P(sost,1) + P(sost,2))
        Res(n,i) = 3;
        sost = 3;
  end
  end
  end
end
end

disp('');
disp('Результат: ');
disp(Res);

disp('------------------------------------------');
disp('Задача 3. Моделирование времени блуждания точки в области.');
disp(' ');

p = [0.15 0.15 0.35 0.35]
A = [2 6]
B = [-2 6]
C = [-2 -4]
D = [2 -4]

n = 100
allTime = 0;
for i=1:n
    coord = [0 0];
    k = 0;
    while(((-2<=coord(1))&&(coord(1)<=2))&&((-4<=coord(2))&&(coord(2)<=6)))
        alpha3(i) = rand(1);
        if(alpha3(i)<=p(1))
            coord(1) = coord(1) + 1;
        else if((p(1)<alpha3(i))&&(alpha3(i)<=p(1)+p(2)))
            coord(2) = coord(2) + 1;
        else if((p(1)+p(2)<alpha3(i))&&(alpha3(i)<=p(1)+p(2)+p(3)))
             coord(1) = coord(1) - 1;   
        else if(alpha3(i)>p(1)+p(2)+p(3))
             coord(2) = coord(2) - 1;       
        end
        end
        end
        end
    k = k + 1;
    end
    allTime = allTime + k-1;
    disp(sprintf('%g  ;  (%g,%g)', k-1,coord(1),coord(2)));

end
disp('Среднее время:')
disp(allTime / n)


    

In [ ]: