И математическая статистика
ГБОУ Колледж «ЦАРИЦЫНО» Отделение управления и информационных технологий
Теория Вероятностей И математическая статистика
Тема 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) )) Результат.
Популярное: ©2015-2024 megaobuchalka.ru Все материалы представленные на сайте исключительно с целью ознакомления читателями и не преследуют коммерческих целей или нарушение авторских прав. (597)
|
Почему 1285321 студент выбрали МегаОбучалку... Система поиска информации Мобильная версия сайта Удобная навигация Нет шокирующей рекламы |