taudata Analytics
Random Number Generator (RNG) & Simulasi (Monte Carlo)
https://taudata.blogspot.com/2022/04/PyRNG.html
Video Penjelasan: https://www.youtube.com/c/taudataAnalytics
Outline RNG: ¶
- Introduction to Random Number [True Randomness]
- Applications, Libraries, and Programs
- PseudoRandom
- Statistically Distributed Random Numbers Generator A. Invers Fungsi Distribusi B. Metode Penolakan
Referensi:
# Importing Modules untuk Notebook ini
import warnings; warnings.simplefilter('ignore')
import matplotlib.pyplot as plt, seaborn as sns, numpy as np
sns.set(style="ticks", color_codes=True)
"Done"
'Done'
Pendahuluan – True Randomness ¶
Dari definisi Sebenarnya True Random Number Generator (TRNG) adalah bilangan acak yang tidak bersifat algoritmik/memiliki aturan tertentu. Sehingga data yang di-generate/hasilkan tidak saling berkorelasi dan bergantung terhadap faktor tertentu (misal waktu, state, seed/nilai awal, dll).
Hal ini hampir tidak mungkin (minimal Sangat Sulit) untuk ditemukan/lakukan dan tidak "praktis" untuk digunakan dalam berbagai aplikasi dalam keseharaian.
Contoh TRNG biasanya pada lab fisika dengan menggunakan special device berdasarkan physical randomness. Namun beberapa tahun ini mulai jamak alat-alat TRNG yang berifat lebih portable (usb devices).
True Random Number Generator - Chaos Example ¶
Pseudo Random Number Generator ¶
- Namun untuk mengenerate bilangan random dengan cara ini sangatlah tidak praktis. Pencatatan setiap even ke data sangatlah tidak efisien dikarenakan simulasi yang menggunakan bilangan random ini biasanya membutuhkan data dalam jumlah yang sangat besar.
- Solusi praktis untuk masalah ini adalah membuat Random Number Generator secara komputasi yang biasa disebut sebagai Pseudo Random number Generator.
Properties of Random Number Generator ¶
Beberapa sifat yang diinginkan dari RNG adalah:
- Randomness: Generator harus lolos uji/tes kerandoman/keacakan yang ketat.
- Long period: Dapat mengenerate data dalam jumlah besar sebelum polanya akan berulang kembali.
- Computational efficiency: Dikarenakan sifat dari masalahnya yang besar, efisiensi komputasi adalah sebuah keharusan.
- Repeatability: Kondisi awal (seed Values) haruslah menentukan sepenuhnya hasil barisan acak/random yang dihasilkan.
- Portability: Barisan bilangan random yang ‘identik’ haruslah dapat dihasilkan dari berbagai komputer yang berbeda dengan seed yang sama.
- Homogeneity: Semua subset dari barisan bilangan random yang dihasilkan haruslah juga random.
Aplikasi RNG ¶
- Simulasi
- Simulasi Model Antrian.
- Simulasi Termodinamik
- Optimasi dari fungsi non-smooth (non differentiable),penyelesaian numerik biasanya dilakukan dengan metode monte carlo (simulated anealing, genetic algorithm).
- Simulasi system dengan resiko tinggi Contoh: Simulasi pada Reaktor Nuklir.
- (Computer/Digital) Security: Encryption, Cryptography.
- Mathematical/Statistical analysis: Stochastic, Monte Carlo, dsb.
- Data Science/Machine Learning.
Libraries - Algorithms ¶
- Rand() LCG adalah library untuk C. Lambat dan tidak mendukung vectorized/concurency.
- SPRNG Didesign khusus untuk parallel computation
- PRNG/Unuran From Uni Salzburg. Widely used, user friendly text based program.
- GNUSL Similar dengan PRNG, tapi cenderung lebih lambat.
- Mersenne Twister High Performance RNG
- Python:
- Python standard Random Library
- Numpy Random number generator
PseudoRandom ¶
- Beberapa ahli komputasi mendefinisikan PSEUDORANDOM Sebagai angka-angka yang di-generate secara algoritmik yang “nampak” independent dan uniform.
- Definisi yang lebih resmi dari pseudorandom adalah barisan angka yang di generate menurut aturan tertentu yang cukup meyakinkan sehingga tidak terdapat tes statistik yang mampu mendeteksi adanya departure (asal) dari kerandomannya.
- Keuntungan dari aturan ini kita mampu memproduksi ulang bilangan-bilangan tersebut untuk tujuan computational cek/reproduction pada komputer lain.
- Dalam pembahasan bilangan acak pseudorandom selanjutnya, maka diasumsikan bilangan yg di hasilkan ($\xi_i$) mempunyai distribusi Standardized rectangular distribution (Uniform):
Teknik Dasar Pembangkitan Bilangan Acak 01 ¶
- Teknik dasar pembangkitan bilangan acak pseudorandom dilakukan dengan menghitung $\xi_i$ dari barisan bilangan integer positive $x_i$.
- Teknik awal yang digunakan untuk membangkitkan bilangan pseudorandom adalah metode midsquare (Metropolis, Von Neumann). Teknik ini menggunakan nilai tengah dari $X_i^2$ sebagai nilai untuk suku berikutnya $(X_{i+1})$.
- Tidak lama berselang (Hull dan Dobell) menemukan bahwa metode midsquare tidak memuaskan, yang kemudian diperbaiki dengan metode Congruential (Lehmer 1951).
- Namun sebelum metode-metode pembangkitan pseudorandom dibahas, perlu diingat bahwa metode-metode ini sebenarnya “lack of theoritical support”, namun tidak ada practitioner yang keberatan dalam menggunakannya selama data yang dihasilkan memenuhi beberapa syarat statistik tertentu.
Teknik Dasar Pembangkitan Bilangan Acak 02 ¶
- Sebuah barisan bilangan pseudorandom dapat dibangkitkan melalui sebuah relasi rekursif: $$ x_i = (ax_{i-1}) \% m $$
- Yaitu suku berikutnya dari sebuah barisan bilangan pseudorandom adalah sisa dari pembagian suku sebelumnya dengan $m$.
- Formulasi ini mempunyai beberapa kelemahan, beberapa diantaranya adalah:
- Barisan bilangannya berulang dengan periode yang selalu $< m$
- Beberapa bilangan mempunyai probabilitas untuk muncul lebih besar
- Formula pseudorandom yang terakhir di generalisasi lebih lanjut oleh Greenberger (1961) dengan: $$ x_i = (ax_{i-1}+c) \% m$$
- Dimana $m$ adalah bilangan yang cukup besar yang ditentukan oleh komputer yang digunakan (biasanya bilangan kelipatan 2 atau 10) dan $a$, $c$, dan $x_i$ adalah bilangan bulat non negative antara 0 dan $m-1$.
- Kedua formulasi terakhir ini disebut metode Congruential, lebih khusus lagi formula Greenberger disebut Multiplicative congruential.
Contoh Perhitungan ¶
$$ x_i = (ax_{i-1}+c)\% m$$- Metode Congruential akan berulang setelah maksimum suku ke m sehingga akan bersifat periodik.
Contoh: $m=16$, $a=3$, $c=1$, dan $x_0=2$.
- Pseudorandom yang dihasilkan oleh formula Greenberger adalah:
- 2, 7, 6, 3, 10, 15, 14, 11, 2, 7, 6, 3, ….
Oleh karena itu kita harus memperhatikan bahwa periodenya adalah lebih besar dari banyaknya data yang hendak kita bangkitkan (generate).
def Greenberger(xo, m, a, c):
return (a*xo + c)%m
m = 16
a = 3
c = 1
xo = 2
x1 = Greenberger(xo,m, a, c)
x1
7
print(xo, end=', ')
x = 2
m = 99999
for i in range(300):
x = Greenberger(x,m, a, c)
print(x, end = ', ')
2, 7, 22, 67, 202, 607, 1822, 5467, 16402, 49207, 47623, 42871, 28615, 85846, 57541, 72625, 17878, 53635, 60907, 82723, 48172, 44518, 33556, 670, 2011, 6034, 18103, 54310, 62932, 88798, 66397, 99193, 97582, 92749, 78250, 34753, 4261, 12784, 38353, 15061, 45184, 35554, 6664, 19993, 59980, 79942, 39829, 19489, 58468, 75406, 26221, 78664, 35995, 7987, 23962, 71887, 15664, 46993, 40981, 22945, 68836, 6511, 19534, 58603, 75811, 27436, 82309, 46930, 40792, 22378, 67135, 1408, 4225, 12676, 38029, 14089, 42268, 26806, 80419, 41260, 23782, 71347, 14044, 42133, 26401, 79204, 37615, 12847, 38542, 15628, 46885, 40657, 21973, 65920, 97762, 93289, 79870, 39613, 18841, 56524, 69574, 8725, 26176, 78529, 35590, 6772, 20317, 60952, 82858, 48577, 45733, 37201, 11605, 34816, 4450, 13351, 40054, 20164, 60493, 81481, 44446, 33340, 22, 67, 202, 607, 1822, 5467, 16402, 49207, 47623, 42871, 28615, 85846, 57541, 72625, 17878, 53635, 60907, 82723, 48172, 44518, 33556, 670, 2011, 6034, 18103, 54310, 62932, 88798, 66397, 99193, 97582, 92749, 78250, 34753, 4261, 12784, 38353, 15061, 45184, 35554, 6664, 19993, 59980, 79942, 39829, 19489, 58468, 75406, 26221, 78664, 35995, 7987, 23962, 71887, 15664, 46993, 40981, 22945, 68836, 6511, 19534, 58603, 75811, 27436, 82309, 46930, 40792, 22378, 67135, 1408, 4225, 12676, 38029, 14089, 42268, 26806, 80419, 41260, 23782, 71347, 14044, 42133, 26401, 79204, 37615, 12847, 38542, 15628, 46885, 40657, 21973, 65920, 97762, 93289, 79870, 39613, 18841, 56524, 69574, 8725, 26176, 78529, 35590, 6772, 20317, 60952, 82858, 48577, 45733, 37201, 11605, 34816, 4450, 13351, 40054, 20164, 60493, 81481, 44446, 33340, 22, 67, 202, 607, 1822, 5467, 16402, 49207, 47623, 42871, 28615, 85846, 57541, 72625, 17878, 53635, 60907, 82723, 48172, 44518, 33556, 670, 2011, 6034, 18103, 54310, 62932, 88798, 66397, 99193, 97582, 92749, 78250, 34753, 4261, 12784, 38353, 15061, 45184, 35554, 6664, 19993, 59980, 79942, 39829, 19489, 58468, 75406, 26221, 78664, 35995, 7987, 23962, 71887, 15664, 46993, 40981, 22945, 68836,
Bilangan Random Uniform (0,1) Menggunakan Formula Dasar ¶
# Sengaja dibuat looping (bukan vectorized) untuk memudahkan dalam mengerti prosesnya.
# Silahkan latihan dengan merubahnya dalam bentuk vectorized.
def uniform(mu=16807, mod_=(2**31)-1, seed=11, size=1):
U = np.zeros(size)
x = (seed*mu+1)%mod_
U[0] = x/mod_
for i in range(1,size):
x = (x*mu+1)%mod_
U[i] = x/mod_
return U
uniform(size=20)
array([8.60905275e-05, 4.46923496e-01, 4.43204159e-01, 9.32302388e-01, 2.06240681e-01, 2.87119356e-01, 6.15018271e-01, 6.12083995e-01, 2.95697290e-01, 7.84355087e-01, 6.55948980e-01, 5.34508298e-01, 4.80963384e-01, 5.51603139e-01, 7.93960019e-01, 8.60457956e-02, 1.71687087e-01, 5.44873295e-01, 6.85463591e-01, 5.86571696e-01])
Visualisasi ¶
X = uniform(size=10**2)
plot = sns.displot(x=X, kde=True)
X = uniform(size=10**5)
plot = sns.displot(x=X, kde=True)
Bilangan Random Uniform (0,1) Menggunakan Module Standard Python ¶
# seed the pseudorandom number generator
from random import seed
from random import random
seed(1) # seed random number generator
print(random(), random(), random()) # generate some random numbers
seed(7)
print(random(), random(), random())
seed(7)
print(random(), random(), random())
0.13436424411240122 0.8474337369372327 0.763774618976614 0.32383276483316237 0.15084917392450192 0.6509344730398537 0.32383276483316237 0.15084917392450192 0.6509344730398537
# Random Integer?
from random import randint
# generate some integers
for _ in range(10):
value = randint(0, 10)
print(value, end=', ')
1, 8, 1, 5, 9, 0, 8, 3, 0, 1,
# Memilih secara acak dari list
from random import choice
seed(1)
List_ = [1, 6, 3, 88, 45, 67, 21, 90, 93]
print(List_)
# make choices from the sequence
for _ in range(5):
selection = choice(List_)
print(selection, end=', ')
[1, 6, 3, 88, 45, 67, 21, 90, 93] 3, 6, 45, 6, 90,
# Mengacak List
from random import shuffle
seed(1)
List_ = [1, 6, 3, 88, 45, 67, 21, 90, 93]
print(List_)
shuffle(List_) # hati-hati defaultnya inplace Operation
print(List_)
[1, 6, 3, 88, 45, 67, 21, 90, 93] [67, 21, 90, 45, 88, 1, 93, 6, 3]
Bilangan Random Uniform (0,1) Menggunakan Module Numpy ¶
# seed the pseudorandom number generator
from numpy.random import seed
from numpy.random import rand
seed(1)
print(rand(3))
seed(7)
print(rand(3))
seed(7)
print(rand(3))
[4.17022005e-01 7.20324493e-01 1.14374817e-04] [0.07630829 0.77991879 0.43840923] [0.07630829 0.77991879 0.43840923]
from numpy.random import randint
seed(1)
values = randint(0, 10, 20)
print(values)
[5 8 9 5 0 0 1 7 6 9 2 4 5 2 4 2 4 7 7 9]
from numpy.random import shuffle
seed(1)
List_ = [1, 6, 3, 88, 45, 67, 21, 90, 93]
print(List_)
shuffle(List_)
print(List_)
[1, 6, 3, 88, 45, 67, 21, 90, 93] [93, 3, 21, 90, 6, 1, 45, 88, 67]
Statistically Distributed RNG ¶
- Biasanya suatu simulasi kurang representative karena data aslinya non stationer (bergantung waktu atau urutan) dan autocorrelate (tidak saling bebas antar observasi), maka variabel random di generate tidak hanya random uniform tapi variable tersebut haruslah memiliki distribusi tertentu.
metode 01: Invers Fungsi Distribusi¶
Syarat: X variabel acak kontinu, strictly increasing di $0\leq F(X)\leq 1$.
- Algoritma:
- Generate U = U(0,1)
- $X = F^{-1}(U)$
- Bukti: $P(X\leq x) = P(F^{-1}(U)\leq x) = P(U\leq F(x)) = F(x)$
Contoh Invers Fungsi Distribusi ¶
Misal $X$ berdistribusi eksponensial dengan mean=$\beta$ maka: $$ F(x)= \left\{ \begin{array}{ll} 1-e^{-x/\beta}, \ \ \ x \geq 0 \\ 0, \ \ \ lainnya \\ \end{array} \right. $$
Sehingga $F^{-1}(x)= -\beta ln(1-u)$ atau $X= –\beta ln(u)$
karena $U(0,1)$ maka $U=1–U$.
a = -2
b = 2
N = 10**1
seed(1)
U = rand(N)
x = [a+(b-a)*u for u in U]
x
[-0.331911981189704, 0.8812979737686324, -1.9995425007306205, -0.7906697094726409, -1.4129764367315478, -1.6306456209248088, -1.2549591544893164, -0.617757091827809, -0.41293010307732025, 0.15526693601342778]
# Visualisasi
N = 10**5
seed(1)
U = rand(N)
x = [a+(b-a)*u for u in U]
plot = sns.displot(x=x, kde=True)
beta = 5
N = 10**4
seed(1)
U = rand(N)
x = [-beta*np.log(u) for u in U]
x[:10]
[4.3730814475232656, 1.6402674158191894, 45.38014815677231, 5.981138168201195, 9.594923396126656, 11.911465400154045, 8.403052985182658, 5.312934427688155, 4.62202438557981, 3.0918988850776294]
# Visualisasi
plot = sns.displot(x=x, kde=True)
N = 10**4
seed(1)
def randNormal(N):
result = []
while len(result)<N:
U1, U2 = rand(), rand()
V1, V2 = 2*U1-1, 2*U2-1
W = V1**2 + V2**2
if W<=1:
y = np.sqrt((-2*np.log(W)/W))
x1, x2 = V1*y, V2*y
result.extend([x1, x2])
return np.array(result)
x = randNormal(N)
x[:10]
array([-0.61175641, 1.62434536, -1.07296862, -0.52817175, -2.3015387 , 0.86540763, -0.7612069 , 1.74481176, -0.24937038, 0.3190391 ])
# Visualisasi
plot = sns.displot(x=x, kde=True)
# Random Normal Numpy
x = np.random.randn(N)
x[:10]
array([-0.12247391, 0.22816982, -0.35230513, -0.83055344, -0.26108982, 0.16935423, 0.6736231 , -0.32720161, -0.30529915, 0.52486533])
import seaborn as sns
data = sns.load_dataset('tips')
data = data[['total_bill', 'tip', 'size']].values
data
array([[16.99, 1.01, 2. ], [10.34, 1.66, 3. ], [21.01, 3.5 , 3. ], [23.68, 3.31, 2. ], [24.59, 3.61, 4. ], [25.29, 4.71, 4. ], [ 8.77, 2. , 2. ], [26.88, 3.12, 4. ], [15.04, 1.96, 2. ], [14.78, 3.23, 2. ], [10.27, 1.71, 2. ], [35.26, 5. , 4. ], [15.42, 1.57, 2. ], [18.43, 3. , 4. ], [14.83, 3.02, 2. ], [21.58, 3.92, 2. ], [10.33, 1.67, 3. ], [16.29, 3.71, 3. ], [16.97, 3.5 , 3. ], [20.65, 3.35, 3. ], [17.92, 4.08, 2. ], [20.29, 2.75, 2. ], [15.77, 2.23, 2. ], [39.42, 7.58, 4. ], [19.82, 3.18, 2. ], [17.81, 2.34, 4. ], [13.37, 2. , 2. ], [12.69, 2. , 2. ], [21.7 , 4.3 , 2. ], [19.65, 3. , 2. ], [ 9.55, 1.45, 2. ], [18.35, 2.5 , 4. ], [15.06, 3. , 2. ], [20.69, 2.45, 4. ], [17.78, 3.27, 2. ], [24.06, 3.6 , 3. ], [16.31, 2. , 3. ], [16.93, 3.07, 3. ], [18.69, 2.31, 3. ], [31.27, 5. , 3. ], [16.04, 2.24, 3. ], [17.46, 2.54, 2. ], [13.94, 3.06, 2. ], [ 9.68, 1.32, 2. ], [30.4 , 5.6 , 4. ], [18.29, 3. , 2. ], [22.23, 5. , 2. ], [32.4 , 6. , 4. ], [28.55, 2.05, 3. ], [18.04, 3. , 2. ], [12.54, 2.5 , 2. ], [10.29, 2.6 , 2. ], [34.81, 5.2 , 4. ], [ 9.94, 1.56, 2. ], [25.56, 4.34, 4. ], [19.49, 3.51, 2. ], [38.01, 3. , 4. ], [26.41, 1.5 , 2. ], [11.24, 1.76, 2. ], [48.27, 6.73, 4. ], [20.29, 3.21, 2. ], [13.81, 2. , 2. ], [11.02, 1.98, 2. ], [18.29, 3.76, 4. ], [17.59, 2.64, 3. ], [20.08, 3.15, 3. ], [16.45, 2.47, 2. ], [ 3.07, 1. , 1. ], [20.23, 2.01, 2. ], [15.01, 2.09, 2. ], [12.02, 1.97, 2. ], [17.07, 3. , 3. ], [26.86, 3.14, 2. ], [25.28, 5. , 2. ], [14.73, 2.2 , 2. ], [10.51, 1.25, 2. ], [17.92, 3.08, 2. ], [27.2 , 4. , 4. ], [22.76, 3. , 2. ], [17.29, 2.71, 2. ], [19.44, 3. , 2. ], [16.66, 3.4 , 2. ], [10.07, 1.83, 1. ], [32.68, 5. , 2. ], [15.98, 2.03, 2. ], [34.83, 5.17, 4. ], [13.03, 2. , 2. ], [18.28, 4. , 2. ], [24.71, 5.85, 2. ], [21.16, 3. , 2. ], [28.97, 3. , 2. ], [22.49, 3.5 , 2. ], [ 5.75, 1. , 2. ], [16.32, 4.3 , 2. ], [22.75, 3.25, 2. ], [40.17, 4.73, 4. ], [27.28, 4. , 2. ], [12.03, 1.5 , 2. ], [21.01, 3. , 2. ], [12.46, 1.5 , 2. ], [11.35, 2.5 , 2. ], [15.38, 3. , 2. ], [44.3 , 2.5 , 3. ], [22.42, 3.48, 2. ], [20.92, 4.08, 2. ], [15.36, 1.64, 2. ], [20.49, 4.06, 2. ], [25.21, 4.29, 2. ], [18.24, 3.76, 2. ], [14.31, 4. , 2. ], [14. , 3. , 2. ], [ 7.25, 1. , 1. ], [38.07, 4. , 3. ], [23.95, 2.55, 2. ], [25.71, 4. , 3. ], [17.31, 3.5 , 2. ], [29.93, 5.07, 4. ], [10.65, 1.5 , 2. ], [12.43, 1.8 , 2. ], [24.08, 2.92, 4. ], [11.69, 2.31, 2. ], [13.42, 1.68, 2. ], [14.26, 2.5 , 2. ], [15.95, 2. , 2. ], [12.48, 2.52, 2. ], [29.8 , 4.2 , 6. ], [ 8.52, 1.48, 2. ], [14.52, 2. , 2. ], [11.38, 2. , 2. ], [22.82, 2.18, 3. ], [19.08, 1.5 , 2. ], [20.27, 2.83, 2. ], [11.17, 1.5 , 2. ], [12.26, 2. , 2. ], [18.26, 3.25, 2. ], [ 8.51, 1.25, 2. ], [10.33, 2. , 2. ], [14.15, 2. , 2. ], [16. , 2. , 2. ], [13.16, 2.75, 2. ], [17.47, 3.5 , 2. ], [34.3 , 6.7 , 6. ], [41.19, 5. , 5. ], [27.05, 5. , 6. ], [16.43, 2.3 , 2. ], [ 8.35, 1.5 , 2. ], [18.64, 1.36, 3. ], [11.87, 1.63, 2. ], [ 9.78, 1.73, 2. ], [ 7.51, 2. , 2. ], [14.07, 2.5 , 2. ], [13.13, 2. , 2. ], [17.26, 2.74, 3. ], [24.55, 2. , 4. ], [19.77, 2. , 4. ], [29.85, 5.14, 5. ], [48.17, 5. , 6. ], [25. , 3.75, 4. ], [13.39, 2.61, 2. ], [16.49, 2. , 4. ], [21.5 , 3.5 , 4. ], [12.66, 2.5 , 2. ], [16.21, 2. , 3. ], [13.81, 2. , 2. ], [17.51, 3. , 2. ], [24.52, 3.48, 3. ], [20.76, 2.24, 2. ], [31.71, 4.5 , 4. ], [10.59, 1.61, 2. ], [10.63, 2. , 2. ], [50.81, 10. , 3. ], [15.81, 3.16, 2. ], [ 7.25, 5.15, 2. ], [31.85, 3.18, 2. ], [16.82, 4. , 2. ], [32.9 , 3.11, 2. ], [17.89, 2. , 2. ], [14.48, 2. , 2. ], [ 9.6 , 4. , 2. ], [34.63, 3.55, 2. ], [34.65, 3.68, 4. ], [23.33, 5.65, 2. ], [45.35, 3.5 , 3. ], [23.17, 6.5 , 4. ], [40.55, 3. , 2. ], [20.69, 5. , 5. ], [20.9 , 3.5 , 3. ], [30.46, 2. , 5. ], [18.15, 3.5 , 3. ], [23.1 , 4. , 3. ], [15.69, 1.5 , 2. ], [19.81, 4.19, 2. ], [28.44, 2.56, 2. ], [15.48, 2.02, 2. ], [16.58, 4. , 2. ], [ 7.56, 1.44, 2. ], [10.34, 2. , 2. ], [43.11, 5. , 4. ], [13. , 2. , 2. ], [13.51, 2. , 2. ], [18.71, 4. , 3. ], [12.74, 2.01, 2. ], [13. , 2. , 2. ], [16.4 , 2.5 , 2. ], [20.53, 4. , 4. ], [16.47, 3.23, 3. ], [26.59, 3.41, 3. ], [38.73, 3. , 4. ], [24.27, 2.03, 2. ], [12.76, 2.23, 2. ], [30.06, 2. , 3. ], [25.89, 5.16, 4. ], [48.33, 9. , 4. ], [13.27, 2.5 , 2. ], [28.17, 6.5 , 3. ], [12.9 , 1.1 , 2. ], [28.15, 3. , 5. ], [11.59, 1.5 , 2. ], [ 7.74, 1.44, 2. ], [30.14, 3.09, 4. ], [12.16, 2.2 , 2. ], [13.42, 3.48, 2. ], [ 8.58, 1.92, 1. ], [15.98, 3. , 3. ], [13.42, 1.58, 2. ], [16.27, 2.5 , 2. ], [10.09, 2. , 2. ], [20.45, 3. , 4. ], [13.28, 2.72, 2. ], [22.12, 2.88, 2. ], [24.01, 2. , 4. ], [15.69, 3. , 3. ], [11.61, 3.39, 2. ], [10.77, 1.47, 2. ], [15.53, 3. , 2. ], [10.07, 1.25, 2. ], [12.6 , 1. , 2. ], [32.83, 1.17, 2. ], [35.83, 4.67, 3. ], [29.03, 5.92, 3. ], [27.18, 2. , 2. ], [22.67, 2. , 2. ], [17.82, 1.75, 2. ], [18.78, 3. , 2. ]])
# Visualisasi
plot = sns.displot(x=x, kde=True)
Metode Invers Fungsi Distribusi pada Distribusi probabilitas Diskrit ¶
- Generate U(0,1)
- Cari i terkecil sehingga u≤F(xi), X= xi
Contoh¶
- misal variabel random 1,2,3,4 akan digenerate dengan pmf: 1/6, 1/3, 1/3, dan 1/6
- Generate U
$$ F(U)= \left\{ \begin{array}{ll} U < 1/6, \ \ \ maka \ \ x=1\\ 1/6\leq U< 1/2, \ \ \ maka \ \ x=2\\ 1/2\leq U< 5/6, \ \ \ maka \ \ x=3\\ 5/6 \leq U, \ \ \ maka \ \ x=4\\ \end{array} \right. $$
N = 10**4
seed(1)
U = rand(N)
for i,u in enumerate(U):
if u<1/6:
U[i] = 1
elif u>=1/6 and u<0.5:
U[i] = 2
elif u>=1/2 and u<5/6:
U[i] = 3
else:
U[i] = 4
X = np.array(U)
X[:10]
array([2., 3., 1., 2., 1., 1., 2., 2., 2., 3.])
plt.figure(figsize=(8,6))
p = sns.countplot(x=X)
N = 10**4
p = 0.5
seed(1)
U = rand(N)
for i,u in enumerate(U):
if u<=p:
U[i] = 1
else:
U[i] = 0
X = np.array(U)
X[:10]
array([1., 0., 1., 1., 1., 1., 1., 1., 1., 0.])
plt.figure(figsize=(8,6))
p = sns.countplot(x=X)
from tqdm import tqdm
N = 10**4
p = 0.5
t = 10**3
def bernoulli(p, t):
U = rand(t)
for i,u in enumerate(U):
if u<=p:
U[i] = 1
else:
U[i] = 0
return U
seed(1)
X = np.zeros(N)
for i,x in tqdm(enumerate(X)):
X[i] = np.sum(bernoulli(p, t))
X[:10]
10000it [00:03, 2998.30it/s]
array([494., 471., 523., 486., 511., 461., 520., 508., 523., 493.])
plt.figure(figsize=(8,6))
p = sns.countplot(x=X)
Monte Carlo ¶
- Metode Monte Carlo adalah sebuah teknik yang menggunakan random numbers dan probabilitas untuk menyelesaikan suatu masalah.
- Kata Monte Carlo pertama kali di pelopori oleh S. Ulam dan Nicholas Metropolis berdasarkan permainan poker yang terkenal di Monte Carlo, Monaco (Hoffman, 1998; Metropolis and Ulam, 1949).
- Metode Monte Carlo biasanya baik untuk digunakan ketika sebuah probabilitas kemunculan suatu proses jelas sifatnya namun hasil outputnya sulit untuk ditentukan (Poisson, Binomial, Normal, dll). Monte Carlo biasa digunakan pada masalah yang secara analitis sulit atau tidak dapat untuk diselesaikan.
- Secara umum model Matematika dibagi menjadi dua macam, yaitu deterministik dan stochastic.
Deterministik¶
- Simulasi Komputer menggunakan model komputasi untuk meniru (imitate) masalah pada aplikasi dunia nyata atau membuat prediksi. Ketika kita membuat model diperlukan beberapa parameter input untuk modelnya dan beberapa persamaan yang menggunakan input-input tersebut untuk menghasilkan himpunan output (response Variables). Model seperti ini biasanya disebut deterministik, yaitu kita pasti akan mendapatkan hasil yang sama betapapun kita mengulang perhitungannya. Input tertentu menghasilkan output yang tertentu pula.
Contoh Model Deterministik ¶
Simulasi Monte Carlo - 01 ¶
- Adalah sebuah metode yang secara iteratif menghitung model deterministik menggunakan bilangan random sebagai inputnya. Metode ini biasanya digunakan ketika modelnya komplex, nonlinear, atau menyangkut beberapa parameter yang tidak pasti (uncertain). Sebuah simulasi monte carlo biasanya dilakukan lebih dari setidaknya ribuan kali dalam menghitung modelnya, sehingga hardware dan teknik pemrogramannya haruslah mendukung hal tersebut.
- Dengan menggunakan input yang random, sebenarnya kita telah merubah model dari deterministik ke stochastic.
- Metode Monte Carlo adalah salah satu metode yang menganalisa suatu perkembangan dari ketidakpastian (uncertainty propagation), dimana tujuan utamanya adalah menentukan seberapa random variasinya (random variation), kurangnya pengetahuan (lack of knowledge), atau error yang mempengaruhi sensitivitas, performa, atau reliabilitas dari sistem yang dimodelkan.
Simulasi Monte Carlo -02 ¶
Metode Monte Carlo dikategorikan sebagai metode Sampling karena inputnya di generate secara random dari suatu distribusi kepekatan tertentu (PDF) untuk mensimulasikan suatu teknik sampling dari populasi yang sesungguhnya. Sehingga kita berusaha untuk memilih distribusi inputnya yang sedekat mungkin dengan data yang telah dimiliki atau paling tidak yang paling representatif dari pengetahuan kita saat ini.
Data yang dihasilkan dari simulasinya dapat ditampilkan sebagai distribusi probabilitas (atau Histogram), atau di konversi ke bar error, reliabilitas prediksi, zone toleransi, dan Interval konfidensi.
5 LANGKAH SEDERHANA METODE MONTE CARLO ¶
- Buat Model Parametriknya, y = f(x1, x2, ..., xq).
- Generate Input random, xi1, xi2, ..., xiq.
- Hitung model dan simpan hasilnya sebagai yi.
- Ulangi steps 2 and 3 for i = 1 to n.
- Analisa hasil dengan histograms, summary statistics, confidence intervals, etc.
CONTOH PERMASALAHAN SEDERHANA YANG DIPECAHKAN DENGAN METODE MONTE CARLO ¶
- Misal terdapat sebuah bujur sangkar dengan panjang sisinya 1. Di dalam bujur sangkar tsb terdapat ¼ lingkaran dengan jari-jari 1, sehingga luasnya adalah pi/4. Dengan demikian sebuah titik (x,y) didalam bujur sangkar tetapi diluar lingkaran akan memenuhi pertidaksamaan x2+y2>1.
- Dengan Monte carlo kita akan coba aproksimasi nilai pi dengan sistem di atas.
- Kita akan men-generate N bilangan Random uniform [0,1] = (x,y) dan memeriksa apakah titik tersebut berada di dalam lingkaran. Ratio Jumlah titik di dalam lingkaran dan total titik yang digunakan akan mendekati ratio luas lingkaran dan bujur sangkar. Sehingga nilai pi dapat diaproksimasi dengan formula:
- 4*(jumlah titik dlm Ligkaran)/(Total titik.)
# Code sengaja tidak dibuat efisien (vectorized) agar lebih mudah dipahami
Nrand = 10*6
NInside = 0
seed(1)
for i in range(Nrand):
if rand()**2 + rand()**2 <= 1.0:
NInside += 1
piapprox = 4*NInside/Nrand
print("pi asli = {}, pendekatan MC={}, error={}".format(np.pi, piapprox, abs(piapprox-np.pi)))
pi asli = 3.141592653589793, pendekatan MC=3.066666666666667, error=0.07492598692312624
Hit Or Miss MONTE CARLO ¶
Latihan Hit Or Miss MONTE CARLO ¶
# Jawaban Latihan 1, integral dari 0 ke 1 dari f(x) thd x dengan f(x)=[exp(x)-1]/[exp(1) - 1]
def f(x):
return (np.exp(x)-1)/(np.exp(1)-1)
seed(1)
N=50
teta = 0.418
sigma = 0.286 #Solusi Eksak teta dan sigma
xrand = rand(N)
yrand = rand(N) # Random Uniform x dan y
print('$$$ Hit or Miss Monte Carlo $$$') # Hit or Miss Monte Carlo
fxrand=f(xrand)
CheckValue = fxrand>=yrand
Hit_or_Miss=sum(CheckValue)/N
print('Hasil Approksimasi = ', Hit_or_Miss)
print('Error = ', abs(teta-Hit_or_Miss))
print('Standar Error = ', np.sqrt(teta*(1-teta)/N))
print('Aproksimasi Error = ', np.sqrt(Hit_or_Miss*(1-Hit_or_Miss)/N))
print('Confidence Interval 95% = ', Hit_or_Miss, ' +/- ', 2*np.sqrt(Hit_or_Miss*(1-Hit_or_Miss)/N))
$$$ Hit or Miss Monte Carlo $$$ Hasil Approksimasi = 0.42 Error = 0.0020000000000000018 Standar Error = 0.06975327949279518 Aproksimasi Error = 0.06979971346646059 Confidence Interval 95% = 0.42 +/- 0.13959942693292118
Tidak ada komentar:
Posting Komentar
Relevant & Respectful Comments Only.