In [86]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import random
import sympy as sp
import pandas as pd
sp.init_printing(use_unicode=True, wrap_line=False, no_global=True)
np.set_printoptions(precision=3)
rnd = random.Random('123')
rand = rnd.random
In [2]:
print('Задача 1. Моделирование непрерывной случайной величины с помощью обратной функции.')
print('Дана плотность распределения')
print(' ')
# Строим график плотности распределения
# Количество реализаций
num = 50
a = 2.0
b = 5.0
# Задаём интервала по x
x1 = np.linspace(-10, 3)
x1_2 = np.zeros(num)
x1_2[:] = 2.0
x2 = np.linspace(a, b)
x2_3 = np.zeros(num)
x2_3[:] = 5.0
x3 = np.linspace(5, 15)
x1_2
# Считаем значения функций на всех интервалах
y1 = np.zeros(num)
y1[:] = 0.0
y1_2 = np.linspace(0, (1 / 39) * (x1_2[0]**2))#0:step_y:(1/63)*x1_2.^2
y2 = np.zeros(num)
y2[:] = (1/39) * (x2**2)
y2_3 = np.linspace(0, (1 / 39) * (x2_3[0]**2)) #0:step_y:(1/63)*x2_3.^2
y3 = np.zeros(num)
#
plt.figure(1)
plt.plot(x1, y1,'r-')
plt.plot(x1_2, y1_2,'r-')
plt.plot(x2, y2,'r-')
plt.plot(x2_3, y2_3,'r-')
plt.plot(x3, y3,'r-')
plt.ylabel('f_ksi(x)')
plt.xlabel('x')
plt.hold(True)
plt.grid(True)
x = sp.Symbol('x')
F_x = sp.integrate((1/39) * (x**2), (x, 3, x))
print('F_x:')
sp.pretty_print(F_x)
In [7]:
print(sp.cancel(F_x))
arg_x = 0.00529100529100529
arg_b = 0.142857142857143
In [3]:
# Строим график функции распределения
# Задаём интервала по x
x1 = np.linspace(-10, 2)
x2 = np.linspace(2, 5)
x3 = np.linspace(5, 15)
# Считаем значения функций на всех интервалах
y1 = np.zeros(x1.shape)
y2 = np.zeros(x2.shape)
y2[:] = (1/117) * (x2**3-8)
y3 = np.ones(x3.shape)
plt.figure(1)
plt.plot(x1, y1,'r-')
plt.plot(x2, y2,'r-')
plt.plot(x3, y3,'r-')
plt.ylabel('F_ksi(x)')
plt.xlabel('x')
plt.grid(True)
plt.hold(True)
print('-------------------------- ')
# Количество реализаций
print('Количество реализаций n:')
n1 = 100
print(n1)
# # Массив псевдослучайных чисел
print('Массив псевдослучайных чисел альфа:')
genalpha1 = np.random.uniform(size=100) # rand(100,1)
print(genalpha1)
# Получим n-реализаций случайной величины кси с помощью моделирующей
# Массив n-реализаций случайной величины кси A
print('Массив n-реализаций случайной величины кси A:')
A1 = ((117) * genalpha1 + 8) ** (1/3)
print(A1)
In [19]:
print('-------------------------- ')
print('Задача 2. Моделирование непрерывной случайной величины с кусочно-линейной плотностью.')
print('Дана плотность распределения.')
print(' ')
A = np.array([-3, 0.05])
B = np.array([-1, 0.05])
C = np.array([ 2, 0.25])
D = np.array([ 4, 0.05])
E = np.array([7, 0.05])
figure = np.array([A, B, C, D, E])
print(figure.shape, figure)
plt.figure(3)
plt.plot(figure[:,0], figure[:, 1], 'r-')
plt.grid(True)
plt.axis('on')
Out[19]:
In [54]:
print('Зададим функцию распределения')
def rand_alpha(a):
if a <= 0.1:
return (a - 0.05) / 0.05
elif 0.1 < a <= 0.55:
return -1.75 + 5 * np.sqrt(-0.0975 + 1.2*a)
elif 0.55 < a <= 0.85:
return 4.5 - 5 * np.sqrt(0.69 - 0.8*a)
else:
return (a - 0.65) / 0.05
print('Создаём 100 реализаций alpha')
alpha = np.random.uniform(size=100)
print(alpha)
print('Для зажаный alpha получаем ksi')
ksi = np.vectorize(rand_alpha)(alpha)
print(ksi)
print('Ksi max:', np.max(ksi))
print('Ksi min:', np.min(ksi))
In [62]:
plt.plot(figure[:,0], figure[:, 1] * 100 - 5, 'r-')
plt.hist(ksi, bins=10)
Out[62]:
In [ ]:
In [67]:
print('---------------------------------')
print('Задача 3. Моделирование трехмерного нормального вектора.')
print(' ')
# Количество реализаций
print('Количество реализаций n:')
n2 = 100
print(n2)
print('Вектор математических ожиданий a.')
a = np.array([-1, 4, 3])
print(a)
print(' ')
print('Ковариационная матрица C.')
C = np.array([
np.array([6, 2, 1]),
np.array([-2, 6, 1]),
np.array([1, 1, 5])
])
print(C)
print(' ')
# создадим матрицу, заполненую нулями.
B = np.zeros((3,3))
B[0,0] = np.sqrt(C[0,0])
print('B11 = ', B[0,0])
B[0,1] = C[0,1] / B[0,0]
print('B12 = ', B[0,1])
B[1,1] = np.sqrt(C[1,1] - B[0,1]**2)
print('B22 = ', B[1,1])
B[0,2] = C[0,2] / B[0,0]
print('B13 = ', B[0,2])
B[1,2] = (C[1,2] - B[0,1] * B[0,2]) / B[1,1]
print('B23 = ', B[1,2])
B[2,2] = np.sqrt(C[2,2] - (B[0,2]**2) - (B[1,2]**2))
print('B33 = ', B[2,2])
print('\r\nВерхнетреугольная матрица B.')
print(B)
In [91]:
print(' ')
print('Найдем ksi.')
# Массив псевдослучайных чисел
# print('Массив псевдослучайных чисел альфа1:')
alpha = np.random.uniform(size=100)
# print(alpha1)
# print('Массив псевдослучайных чисел альфа1:')
alpha2 = np.random.uniform(size=100)
# print(alpha2)
# print('Массив псевдослучайных чисел альфа1:')
alpha3 = np.random.uniform(size=100)
# print(alpha3)
n2 = 100
alpha = np.random.uniform(size=(100,3))
omega = np.zeros((n2, 3))
w = np.zeros((n2, 3))
ksi = np.zeros((n2,3))
for i in range(n2):#i=1:n2
omega[i,0] = 1 - 2*alpha[i, 0]
gamma1 = 1 - 2*alpha[i, 0]
gamma2 = 1 - 2*alpha[i, 1]
D = (gamma1**2) + (gamma2**2)
omega[i,1] = gamma1 * np.sqrt((1-omega[i,0]**2) / D)
omega[i,2] = gamma2 * np.sqrt((1-omega[i,0]**2) / D)
cnt = 0
for j in range(3): #j=1:3:
ksi[i,j] = omega[i,j] * B[j, j] + a[j]
cnt = cnt+1
print('N-реализаций вектора ksi')
# print(ksi)
pksi = pd.DataFrame(ksi, index=None)
pksi.to_excel('lab3_task_3.xlsx')
In [ ]: