И математическая статистика
ГБОУ Колледж «ЦАРИЦЫНО»
Отделение управления и информационных технологий
Теория Вероятностей
И математическая статистика
Тема 1:
Моделирование и визуализация
Случайных данных
Студентка: Зубкова Наталья
Группа: КС 2 – 1
Преподаватель: В.В. Мещеряков
Задание 1.1. Модуль MCG_std.py
#Листинг 1.1. Модуль MCG_std.py
import numpy as np
#x0 - стартовая точка
#N - число генерируемых чисел
def MCG_std(a,c,m,x0,N):
xi = x0
x = np.zeros(N); x = x.astype(np.int32)
r = np.zeros(N)
for k in range(N):
xi = (a*xi + c) % m
x[k] = xi
r[k] = xi/m
return x, r
if __name__ == "__main__":
a = 13; c = 0; m = 31; x0 = 1; N = 12;
#a = 7**5; c = 0; m = 2**31 - 1; x0 = 1; N = 10;
x, r = MCG_std(a,c,m,x0,N)
print('x = ',x)
print('r = ',np.round(r,4))
Результат:
runfile('C:/Users/150114/Desktop/07.12/Модуль MCG_std.py', wdir='C:/Users/150114/Desktop/07.12')
x = [13 14 27 10 6 16 22 7 29 5 3 8]
r = [ 0.4194 0.4516 0.871 0.3226 0.1935 0.5161 0.7097 0.2258 0.9355
0.1613 0.0968 0.2581]
Задание 1.2
#Листинг 1.2. Скрипт MCG_std_subplot.py
import MCG_std as MC
import numpy as np
import matplotlib.pyplot as plt
a = 13; c = 0; m = 31; x0 = 1; N = 100;
#a = 7**5; c = 0; m = 2**31 - 1; x0 = 1; N = 100;
x, r = MC.MCG_std(a,c,m,x0,N);
fig = plt.gcf()
fig.set_size_inches(14,5)
plt.subplot(131)
plt.plot(np.linspace(1,N,N),r, 'b.' )
plt.xlabel('$n$'),
plt.ylabel('$r$')
plt.subplot(132)
r1 = r[0:N/2]; r2 = r[N/2:N + 1]
plt.plot(r1,r2, 'b.' )
plt.xlabel('$r_1$'),
plt.ylabel('$r_2$')
plt.subplot(133)
N = 60
x, r = MC.MCG_std(a,c,m,x0,N);
r1 = r[0:N/2]; r2 = r[N/2:N + 1]
plt.plot(r1,r2, 'b.' )
plt.xlabel('$r_1$'),
plt.ylabel('$r_2$')
Результат:
Мой пример
#Листинг 1.2. Скрипт MCG_std_subplot.py
import MCG_std as MC
import numpy as np
import matplotlib.pyplot as plt
a = 135; c = 8; m = 55; x0 = 1; N = 50;
#a = 7**5; c = 0; m = 2**31 - 1; x0 = 1; N = 100;
x, r = MC.MCG_std(a,c,m,x0,N);
fig = plt.gcf()
fig.set_size_inches(14,5)
plt.subplot(131)
plt.plot(np.linspace(1,N,N),r, 'b.' )
plt.xlabel('$n$'),
plt.ylabel('$r$')
plt.subplot(132)
r1 = r[0:N/2]; r2 = r[N/2:N + 1]
plt.plot(r1,r2, 'b.' )
plt.xlabel('$r_1$'),
plt.ylabel('$r_2$')
plt.subplot(133)
N = 60
x, r = MC.MCG_std(a,c,m,x0,N);
r1 = r[0:N/2];
r2 = r[N/2:N + 1]
plt.plot(r1,r2, 'b.' )
plt.xlabel('$r_1$'),
plt.ylabel('$r_2$')
Результат.
Задача 2.1
#Листинг 2.2. Скрипт MCG_rnd_subplot.py
import MCG_rnd as MC
import numpy as np
import matplotlib.pyplot as plt
low = 10; high = 65; x0 = 6;
fig = plt.gcf()
fig.set_size_inches(14,5)
for n in range(3,6):
plt.subplot(1, 3, n - 2)
N = 10**n
rnd = MC.MCG_rnd(low,high,x0,N)
plt.plot(np.linspace(1,N,N),rnd, 'b.' )
plt.axis([0, N, 10 ,70])
plt.xlabel('$n$'),
plt.ylabel('$rnd$')
Результат.
Задача 3.1
#Листинг 3.1. Скрипт rnd_sum.py
import MCG_rnd as MC
import numpy as np
import matplotlib.pyplot as plt
N = 10**3
low = 0; high = 1;
fig = plt.gcf()
fig.set_size_inches(14,5)
n = np.linspace(1,N,N)
rnd = np.zeros(N)
for k in range(1,5):
plt.subplot(1,4,k)
x0 = 10**k
rnd = rnd + MC.MCG_rnd(low,high,x0,N)
plt.plot(n,rnd, 'b.' )
(plt.text(300,4.3, '$\sum_{k=1}^{k=n}{rn{{d}_{k}}}$, $n =$' + str(k)))
plt.xlabel('$n$')
plt.axis([0, N, -1 ,5])
plt.grid()
Результат.
Задача 4.1
#Листинг 4.1. Скрипт less_than_05_MCG.py
import MCG_rnd as MC
import numpy as np
N = 10**6
print('N =',N)
low = 0; high = 1;
for m in range(1,5):
x0 = 10**(2*m)
rnd = MC.MCG_rnd(low,high,x0,N)
n = 0
for k in range(N):
if rnd[k] < 0.5:
n = n + 1
r = abs(N - 2*n)/N
x0_n_r = np.array([x0,n,r])
print('x0_n_r =',x0_n_r)
Результат.
N = 1000000
x0_n_r = [ 1.00000000e+02 4.98679000e+05 2.64200000e-03]
x0_n_r = [ 1.00000000e+04 5.00736000e+05 1.47200000e-03]
x0_n_r = [ 1.00000000e+06 4.99570000e+05 8.60000000e-04]
x0_n_r = [ 1.00000000e+08 4.99876000e+05 2.48000000e-04]
Задача 5.1
#Листинг 5.1. Скрипт mersenne_random.py
import random
import matplotlib.pyplot as plt
low = 20; high = 80;
fig = plt.gcf()
fig.set_size_inches(14,5)
for k in range(3,6):
plt.subplot(1, 3, k - 2)
N = 10**k
n = range(N)
rnd = [random.uniform(low,high) for i in n]
plt.plot(n, rnd, 'b.')
plt.axis([0, N, 10 ,70])
plt.xlabel( ' $n$ ' ),
plt.ylabel( ' $rnd$ ' )
Результат.
Задача 6.1
6.1. Скрипт less_than_05_Mersenne.py
import numpy as np
N = 10**6
print('N =',N)
low = 0; high = 1;
for m in range(1,5):
x0 = 10**(2*m)
np.random.seed(x0)
rnd = np.random.uniform(low,high,N); #print(rnd)
n = 0
for k in range(N):
if rnd[k] < 0.5:
n = n + 1
r = abs(N - 2*n)/N
x0_n_r = np.array([x0,n,r])
print('x0_n_r =',x0_n_r)
Результат.
N = 1000000
x0_n_r = [ 1.00000000e+02 5.00140000e+05 2.80000000e-04]
x0_n_r = [ 1.00000000e+04 5.00331000e+05 6.62000000e-04]
x0_n_r = [ 1.00000000e+06 4.99994000e+05 1.20000000e-05]
x0_n_r = [ 1.00000000e+08 4.99481000e+05 1.03800000e-03]
Глава 2
Задачи о блуждании по прямой
Задача 7.1
#Листинг 7.1. Модуль steps.py
import numpy as np
def binrnd(n):
rnd = (np.append(np.array([0]),
np.random.random_integers(0, 1, n)))
return rnd
def xm(n):
rnd = binrnd(n);
rnd_temp = np.array(rnd)
x = rnd_temp
for k in range(n):
x[k + 1] = x[k + 1] + x[k]
return x, rnd
if __name__ == "__main__":
n = 10
x, rnd = xm(n)
print('rnd =',rnd)
print('x =',x)
Результат.
rnd = [0 0 0 0 0 0 1 1 1 1 1]
x = [0 0 0 0 0 0 1 2 3 4 5]
Задача 7.2
#Листинг 7.2. Скрипт random_walk.py
import steps
import numpy as np
import matplotlib.pyplot as plt
n = 3; N = 2
# Графическое окно
plt.close('all')
fig = plt.gcf()
fig.set_size_inches(5,5)
plt.xlabel('$t$');
plt.ylabel('$x$')
plt.grid()
for j in range(N):
x, rnd = steps.xm(n)
t = np.linspace(0,n,n + 1);
t = t.astype(np.int32)
plt.plot(t,x,'ro-')
if n <= 20:
plt.xlim(-.5,n + 0.5);
plt.ylim(-.5,n + 0.5)
plt.xticks(np.arange(0,n + 1,1))
plt.yticks(np.arange(0,n + 1,1))
if N <= 10:
print('binrnd',j+1,'=',rnd)
print('x',j+1,'=',x)
Результат.
binrnd 1 = [0 1 0 0]
x 1 = [0 1 1 1]
binrnd 2 = [0 0 0 0]
x 2 = [0 0 0 0]
Задача 8.1
#Листинг 8.1. Скрипт random_walk_subplot.py
import matplotlib.pyplot as plt
import numpy as np
import steps
plt.close('all')
fig = plt.gcf()
fig.set_size_inches(12,4)
N = 5
nmax = 3
for n in range(1,nmax + 1):
for j in range(N):
x, rnd = steps.xm(n)
t = np.linspace(0,n,n + 1); t = t.astype(np.int32)
plt.subplot(1,3,n)
plt.plot(t,x,'ro-')
plt.grid(True)
plt.xlim(-.5,3.5);
plt.ylim(-.5,3.5)
plt.xlabel('$t$');
plt.ylabel('$x$')
plt.xticks(np.arange(0,4,1))
plt.yticks(np.arange(0,4,1))
if n == 1:
plt.text(1.2,0,'1')
plt.text(1.2,1,'1')
elif n == 2:
plt.text(2.2,0,'1')
plt.text(2.2,1,'2')
plt.text(2.2,2,'1')
elif n == 3:
plt.text(3.2,0,'1')
plt.text(3.2,1,'3')
plt.text(3.2,2,'3')
plt.text(3.2,3,'1')
Результат.
#Листинг 8.1. Скрипт random_walk_subplot.py
import matplotlib.pyplot as plt
import numpy as np
import steps
plt.close('all')
fig = plt.gcf()
fig.set_size_inches(12,4)
N = 30
nmax = 3
for n in range(1,nmax + 1):
for j in range(N):
x, rnd = steps.xm(n)
t = np.linspace(0,n,n + 1); t = t.astype(np.int32)
plt.subplot(1,3,n)
plt.plot(t,x,'ro-')
plt.grid(True)
plt.xlim(-.5,3.5);
plt.ylim(-.5,3.5)
plt.xlabel('$t$');
plt.ylabel('$x$')
plt.xticks(np.arange(0,4,1))
plt.yticks(np.arange(0,4,1))
if n == 1:
plt.text(1.2,0,'1')
plt.text(1.2,1,'1')
elif n == 2:
plt.text(2.2,0,'1')
plt.text(2.2,1,'2')
plt.text(2.2,2,'1')
elif n == 3:
plt.text(3.2,0,'1')
plt.text(3.2,1,'3')
plt.text(3.2,2,'3')
plt.text(3.2,3,'1')
Результат.
Задача 8.2
#Листинг 8.2. Скрипт binomial.py
import sympy as sp
a = sp.symbols('a')
b = sp.symbols('b')
s = a + b
for n in range(4):
print(sp.expand(s**n))
результат.
a + b
a**2 + 2*a*b + b**2
a**3 + 3*a**2*b + 3*a*b**2 + b**3
Задача 10.1
#Листинг 10.1. Модуль binopascal.py
import numpy as np
import math
def binopascal(n):
def binomial_coefficient(n, m):
Cnm = (math.factorial(n)/math.factorial(m)/
math.factorial((n - m)))
Cnm = int(Cnm)
return Cnm
def pascal(n):
Cm = np.zeros(n + 1);
Cm = Cm.astype(int);
for m in range(n+1):
Cmn = binomial_coefficient(n, m);
Cm[m] = np.array(Cmn);
return Cm
Cm = pascal(n);
N = int(np.sum(Cm));
Pm = Cm/N;
sumPm = np.sum(Pm);
return Cm, N, Pm, sumPm
if __name__ == "__main__":
n = 3
Cm, N, Pm, sumPm = binopascal(n)
print('Cm =',Cm)
print('N =',N) # N = 2**n, % проверка
print('Pm =',Pm)
print('sumPm =',sumPm)
Результат:
wdir='C:/Users/Acer/Desktop/Моделирование')
Cm = [1 3 3 1]
N = 8
Pm = [ 0.125 0.375 0.375 0.125]
sumPm = 1.0
Задача 10.2
#Листинг 10.2. Скрипт binograph_subplot.py
import binopascal as bin
import matplotlib.pyplot as plt
import numpy as np
fig = plt.gcf()
fig.set_size_inches(10,7)
nmax = 5
for n in range(1,nmax + 1):
Cm, N, Pm, sumPm = bin.binopascal(n);
plt.subplot(nmax/2,2,n)
m = np.arange(n + 1)
width = 0.35
plt.bar(m, Pm, width, align='center', color='red')
plt.axis([-1, 5, 0, 0.6])
plt.xticks(np.arange(0,n + 1,1))
plt.xlabel('$m$')
plt.ylabel('$P_n$')
plt.grid()
Результат:
Глава 3.
Статистический эксперимент Бюффона.
Задача 12.1
#Листинг 12.1. Модуль needle.py
import numpy as np
def needle(l,x0,y0,alpha):
xx1 = x0 - l*np.cos(alpha);
xx2 = x0 + l*np.cos(alpha);
x1 = min(xx1,xx2);
x2 = max(xx1,xx2);
dx = 1e-5;
x = np.arange(x1,x2,dx);
y = y0 + np.tan(alpha)*(x-x0);
return x, y
if __name__ == "__main__":
l,x0,y0,alpha = 0.8, 2.5, 1.5, np.pi/3
x, y = needle(l,x0,y0,alpha)
import numpy as np
import matplotlib.pyplot as plt
plt.axis([0,5,0,3])
plt.plot(x0,y0,'ro')
plt.plot(x,y)
plt.grid()
Результат:
Задача 14.1
#Листинг 14.1. Модуль buffon.py
from numpy.random import uniform
import numpy as np
def buffon(n):
x0 = uniform(0,10,n);
y0 = uniform(0,10,n);
alpha = uniform(0,np.pi,n);
return x0, y0, alpha
if __name__ == "__main__":
n = 5
x0, y0, alpha = buffon(n)
print('x0 =',x0)
print('y0 =',y0)
print('alpha =',alpha)
Результат:
x0 = [ 2.58412323 2.89639669 4.79362979 6.76114968 5.05263329]
y0 = [ 0.46557138 4.93558224 9.66520151 2.72354034 0.22222513]
alpha = [ 2.74291156 2.60389662 0.8735163 1.00904452 2.52620912]
Задача 15.1
Листинг 15.1. Скрипт geometric_probability.py
import sympy as sym
𝜋 = sym.pi
𝛼, l, a = sym.symbols('𝛼 l a')
I = sym.integrate(l*sym.sin(𝛼)/𝜋/a,𝛼)
P2 = I.subs(𝛼,𝜋); P1 = I.subs(𝛼,0)
P = P2-P1
print('P =',P)
Результат:
P = 2*l/(pi*a)
Задача 16.1
#Листинг 16.1. Скрипт buffon_1_run.py
import numpy as np
import matplotlib.pyplot as plt
import buffon, needle, counter
plt.close('all')
fig = plt.gcf()
plt.xlim(-1,11); plt.ylim(-1,11)
fig.set_size_inches(5,5)
plt.xlabel('$x$'); plt.ylabel('$y$')
plt.xticks(np.arange(0,11,2))
plt.yticks(np.arange(0,11,2))
n = 10**0
l = 0.8; a = 1.; xmin = -1; xmax = 11;
Y = np.arange(0.,12.,2*a);
dx = 1e-1; x = np.arange(xmin,xmax,dx);
for k in range(len(Y)):
plt.plot(x,Y[k]*x**0,'r')
x0, y0, alpha = buffon.buffon(n)
for k in range(n):
x, y = needle.needle(l,x0[k],y0[k],alpha[k]);
plt.plot(x,y,'b'),
#print()
#print('Результаты моделирования')
m = counter.counter(n,a,l,y0,alpha); #print('m = ',m)
Pstat = m/n; #print('Pstat = ',Pstat)
Ptheory = 2*l/np.pi/a; #print('Ptheory = ',Ptheory)
Pstat = np.round(Pstat, 4)
Ptheory = np.round(Ptheory, 4)
(plt.title('$n = $' + str(n) + '$, \ m = $' +
str(m) + '$, \ Pstat = $' + str(Pstat) +
'$, \ Ptheory = $' + str(Ptheory) ))
Результат.
