1 Một số khái niệm

Mô phỏng ngẫu nhiên (Random Simulations) hay được gọi với một tên khác

Phương Pháp Monte Carlo (Monte Carlo Method).

Phương pháp mô phỏng Monte Carlo được nhà toán học Stanislaw Ulam

Mô phỏng Là sự bắt chước (hoặc là sự mô tả lại) hoạt động của một đối tượng nào đó mà chúng ta quan tâm hoặc nghiên cứu về nó. Một số ví dụ:

  1. Mô phỏng dòng tiền của một doanh nghiệp;
  2. Mô phỏng kế hoạch marketing.
  2. Mô phỏng quá trình hình thành cục máu đông;
  3. Mô phỏng lượng khí thải của động cơ;
  4. Mô phỏng nợ xấu của một ngân hàng;
  5. Mô phỏng kết quả hoạt động sản xuất kinh doanh của một cửa hàng/Doanh nghiệp; ...

Hệ Thống (system)

Trạng Thái Của Hệ Thống (State of System)

Hệ Thống Rời Rạc (Discrete System)

Hệ Thống Liên Tục (Continuous System)

Hệ Thống Tất Định (Deterministic System)

Hệ Thống Ngẫu Nhiên (Stochastic System)

Mô Hình (Model)

Mô Hình Hoá (Modeling)

2 Các bước cơ bản để thực hiện mô phỏng

  1. Xác định vấn đề cần mô phỏng.
  2. Chọn mô hình và chọn dữ liệu để phục vụ cho việc nghiên cứu.
  3. Xây dựng thuật toán và viết chương trình.
  4. Chạy chương trình.
  5. Kiểm tra.
  6. Phân tích dữ liệu đầu ra.
  7. Kết luận, viết báo cáo.

3 Ưu và nhược điểm của phương pháp mô phỏng

3.1 Benefits of Monte Carlo Simulation

- It enables to make realistic forecasts or manage activities that involve uncertainty.
-It enables to get accurate results by exploring thousands of combinations with “what-if” analysis.
- In this simulation, it’s possible to model interdependent relationships between input variables.
- It helps to improve the quality of decisions.
- It helps to make forecasts for budget, schedule, and other important project work.
- This tool provides graphical and visual data. This helps to improve communication among project team and stakeholders.
- This tool shows the inputs which have the biggest effect on the result.

3.2 Limitations of Monte Carlo Simulation

- It is a time consuming and complicated method.
- It is not an easy method for quantitative risk analysis. You need software to run this method.
- There should be enough samples or inputs to ensure realistic results.
- Results can be effected from the implementors bias.
- Correct analysis and results depend on the quality of the estimates.
- The is tool provides a number of results based on the probabilities. It does not give the actual result.

4 Một Số Ví Dụ Bằng Phần Mềm R

4.1 Tung 1 Con Xúc Xắc

library(tidyverse)
library(DT)
library(ggplot2)
library(ggpattern)
dice <- sample (x =1:6, size = 200,replace = TRUE, prob = c(1/6,1/6,1/6,1/6,1/6,1/6))
table(dice)
## dice
##  1  2  3  4  5  6 
## 40 32 32 30 31 35
dice <- as.data.frame(dice)
ggplot(dice, aes(dice)) +
  geom_bar(fill = 'red') +
  ylab('Số lần xuất hiện') +
  xlab('Số nút')

4.2 Tung 2 Con Xúc Xắc

library(TeachingDemos)
twodice<- dice(200, ndice = 2, plot.it = T)

twodice$sum <- twodice$Red + twodice$Green
ggplot(twodice, aes(sum)) +
  geom_bar(fill = 'red')

4.3 Mô Phỏng Cho Một Số Phân Phối Xác Xuất Thông Dụng

4.3.1 Phân Phối Đều

tmp <- runif(1000, 5, 10)
hist(tmp)

4.3.2 Phân phối Poisson

tmp <- rpois(10, lambda = 1.3)
tmp
##  [1] 2 1 2 2 1 0 1 1 0 0
tmp <- table(tmp)
tmp <- as.data.frame(tmp)
ggplot(tmp, aes(tmp, Freq)) + geom_col(fill = 'red') +
  geom_point() +
  geom_line(aes(as.integer(tmp), Freq), color = 'green')

4.3.3 Phân phối hình học

tmp <- rgeom(10, 0.15)
tmp
##  [1]  5  1  8  4  4  4  2 12 13 10
tmp <- table(tmp)
tmp <- as.data.frame(tmp)
ggplot(tmp, aes(tmp, Freq)) +
  geom_col(fill = 'red') +
  geom_point() +
  geom_line(aes(as.integer(tmp), Freq), color = 'green')

4.3.4 Phân Phối Chuẩn

tmp <- rnorm(100)
hist(tmp)

tmp <- as.data.frame(tmp)
ggplot(data = tmp, aes(tmp)) +
  geom_density(color = 'green')

4.3.5 Phân Phối Mũ

tmp <- rexp(1000,0.5)
hist(tmp)

tmp <- as.data.frame(tmp)
ggplot(data = tmp, aes(tmp)) +
  geom_density(color = 'green')

4.4 Kiểm định một số phân phối thông dụng

Để kiểm định phân phối cho một biến ngẫu nhiên (trong thực tế là 1 biến thống kê) thông thường chúng ta sẽ thực hiện một số bước như sau:

  • Vẽ đồ thị thường là dạng histogram hoặc một số dạng đồ thị theo yêu cầu cụ thể.
  • Dựa vào đồ thị một phần nào đó chúng ta sẽ phát hiện ra phân phối của dãy dữ liệu.
  • Thực hiện các bài toán kiểm định giả thuyết thống kê.
  • Lập lại bước 2 và bước 3 cho đến khi phát hiện ra phân phối của dãy dữ liệu.

4.4.1 Kiểm định phân phối chuẩn

  • Vẽ đồ thị Histogram.
  • Vẽ biểu đồ Q-Q plot.
  • Dùng kiểm định Shapiro-Wilk.
  • Perform a Kolmogorov-Smirnov Test.

Chúng ta sẽ tạo ra 1 biến có phân phối chuẩn và 1 biến có phân phối mũ để đối chứng, sau đó chúng ta sẽ thực hiện lần lượt các thủ tục kiểm định.

nn <- 1000
norr <- rnorm(n = nn, mean = 100, sd = 10)
ex <- 50*rexp(n = nn, rate = 3)
poi <- rpois(n = nn, lambda = 5)

d <- cbind(norr,ex,poi)
d <- as.data.frame(d)
pd <- d |> pivot_longer(cols = c(norr,ex), names_to = 'varName', values_to = 'val')

Vẽ đồ thị histogram

d |> ggplot() +
  geom_histogram(aes(x = norr),binwidth = 1) +
  geom_histogram(aes(x = ex), binwidth = 1, fill = 'red') +
  geom_density(aes(x = norr), color = 'blue')

Vẽ đồ thị QQ plot

pd |> ggplot(aes(sample = val, color = varName)) +
  geom_qq() +
  stat_qq_line()

Chúng ta cũng thực hiện việc vẽ đồ thị như trên nhưng sẽ thực hiện bằng một cách khác.

d |> ggplot() +
  geom_qq(aes(sample = norr), color = 'red') +
  geom_qq_line(aes(sample = norr), color = 'green') +
  geom_qq(aes(sample = ex), color = 'blue') +
  geom_qq_line(aes(sample = ex), color = 'cyan')

Sử dụng kiểm định Shapiro-Wilk Test: Kiểm định này được đề xuất bởi Samuel Sanford Shapiro và Martin Wilk năm 1965.

Với giả thuyết \(H_0: \text{Dãy số liệu có phân phối Chuẩn}\). Để thực hiện thủ tục kiểm định này trong R chúng ta thực hiện như sau:

shapiro.test(d$no)
## 
##  Shapiro-Wilk normality test
## 
## data:  d$no
## W = 0.99837, p-value = 0.475
shapiro.test(d$ex)
## 
##  Shapiro-Wilk normality test
## 
## data:  d$ex
## W = 0.82743, p-value < 2.2e-16
shapiro.test(d$poi)
## 
##  Shapiro-Wilk normality test
## 
## data:  d$poi
## W = 0.96963, p-value = 1.322e-13

Sử dụng kiểm định Kolmogorov-Smirnov (KS): Kiểm định KS là một thủ tục kiểm định phi tham số, kiểm định này được dùng để kiểm định xem một dãy số liệu có một phân phối nào cụ thể nào đó hoặc dùng để kiểm định xem 2 dãy số liệu có được lấy từ cùng một phân phối không?

Để kiểm định xem 1 dãy số liệu có tuân theo một phân phối cụ thể nào không chúng ta thực hiện như sau:

Giả thuyết: \(H_0:\) Dãy số liệu có phân phối cụ thể

ks.test(d$no, y = 'pnorm', mean = 100, sd = 10)
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  d$no
## D = 0.020751, p-value = 0.7823
## alternative hypothesis: two-sided
ks.test(d$ex, y = 'pexp', rate = 3)
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  d$ex
## D = 0.9045, p-value < 2.2e-16
## alternative hypothesis: two-sided
ks.test(d$poi, y = 'pnorm')
## Warning in ks.test.default(d$poi, y = "pnorm"): ties should not be present for
## the Kolmogorov-Smirnov test
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  d$poi
## D = 0.92825, p-value < 2.2e-16
## alternative hypothesis: two-sided

Lưu ý: Tham số thứ 2 y phải là các phân phối liên tục và có thể nhận các giá trị sau: pexp, pnorm.

4.4.2 Goodness of fit test

Đây là một thủ tục kiểm định để kiểm định xem dữ liệu thu thập được có tuân theo một phân phối xác suất nào không. Thủ tục kiểm định này được thực hiện dựa trên ý tưởng: So sánh sự khác biệt giữa giữa dữ liệu quan sát và một phân phối đã biết.

Ví dụ: \(H_0\): there is no significant difference between the observed and the expected frequencies

dice1 <- sample (x =1:6, size = 200,replace = TRUE, prob = c(1/6,1/6,1/6,1/6,1/6,1/6))
dice2 <- sample (x =1:6, size = 200,replace = TRUE, prob = c(1/12,3/12,1/12,1/12,5/12,1/12))
chisq.test(table(dice1), p = c(1/6,1/6,1/6,1/6,1/6,1/6))
## 
##  Chi-squared test for given probabilities
## 
## data:  table(dice1)
## X-squared = 6.1, df = 5, p-value = 0.2966
chisq.test(table(dice1), p = c(1/12,3/12,1/12,1/12,5/12,1/12))
## 
##  Chi-squared test for given probabilities
## 
## data:  table(dice1)
## X-squared = 101.27, df = 5, p-value < 2.2e-16

4.5 Mô phỏng một số quá trình ngẫu nhiên (Stochastic/Random processes)

Quá ngẫu nhiên là một tập hợp các biến ngẫu nhiên được đánh số thứ tự (chỉ số) bằng thời gian. Quá trình ngẫu nhiên được chia làm 2 loại (theo thời gian).

Quá trình ngẫu nhiên rời rạc được ký hiệu \(X = \{X_n, n = 0,1,2,\dots\}\), là một dãy các biến ngẫu nhiên (với mỗi giá trị của \(n, X_n\) là một biến ngẫu nhiên).

Quá trình ngẫu nhiên liên tục được ký hiệu \(X = \{X_t, 0 \le t < \infty\}\), là tập hợp các biến ngẫu nhiên ( với mỗi giá trị của \(t, X_t\) là một biến ngẫu nhiên).

4.5.1 Bước Đi Ngẫu Nhiên (Random Walk)

Dạng đơn giản (Simple random walk):

Cho \(Y_1,Y_2,\dots\) là các biến ngẫu nhiên độc lập cùng phân phối sao cho \(Y_i = 1\) hoặc \(Y_i = -1\). Xét

\[X_k = \sum_{i=1}^kY_i \text{ với }, k \ge 1, X_0 = 0\]

n <- 100
s <- sample(c(1,-1), size = n, prob = c(.5, .5), replace = T)
s[1] <- 0
cs <- numeric(n)
cs <- cumsum(s)
t <- seq(from = 0, to = 99)
d <- as.data.frame(cbind(t,s,cs))
ggplot(d, aes(t,cs)) + geom_point(color = 'blue', size = 1/4) + geom_line(color = 'red')

Dạng Phức Tạp

Cho \(Y_1,Y_2,\dots\) là các biến ngẫu nhiên độc lập cùng phân phối sao cho \(Y_{i+1} = Y_i + X\) với \(X\) là biến ngẫu nhiên có một phân phối nào đó.

\[X_k = \sum_{i=1}^kY_i \text{ với }, k \ge 1, X_0 = 0\]

n <- 100
s <- runif(n, min = -5, max = 5)
s[1] <- 1000
cs <- numeric(n)
cs <- cumsum(s)
t <- seq(from = 0, to = n-1)
d <- as.data.frame(cbind(t,s,cs))
ggplot(d, aes(t,cs)) + geom_point(color = 'blue', size = 1/4) + geom_line(color = 'red')

4.5.2 Quá trình Poisson

Định nghĩa: Cho trước \(\lambda > 0\), Quá trình poisson (hay còn gọi là quá trình đếm) \(\{N(t),t \ge 0\}\) nếu thỏa các yêu điều kiện sau:

  1. \(N(0) = 0\)
  2. Số gia của \(N(t)\) độc lập
  3. Số gia của \(N(t)\) trong một khoảng thời gian \(\tau >0\) có phân phối Poisson với tham số \(\lambda\tau\).

Ví dụ:

t <- 1:50
s <- rpois(length(t) -1, lambda = 8)
s <- c(0,s)
d <- as.data.frame(cbind(t,s))
  
ggplot(d, aes(t,s)) + geom_point(color = 'blue', size = 1/4)+ geom_line(color = 'red')

cs <- cumsum(s)
d <- d |> mutate(cs)
ggplot(d, aes(t,cs)) + geom_point(color = 'blue', size = 1/4) + geom_line(color = 'red')

4.5.3 Chuyển động Brown

Quá trình ngẫu nhiên \(W(t)\), được gọi là chuyển động Brown tiêu chuẩn nếu thỏa các điều kiện sau:

  1. \(W(0) = 0\)
  2. Với mọi \(0 \le t_1 \le t_2\) thì \(W(t_2) - W(t_1) \sim N(0,t_2 - t_1)\).
  3. \(W(t_2) - W(t_1)\) độc lâp.

Ví dụ:

t <- 0:100
sig2 <- 10
x <- rnorm(n = length(t) - 1, sd = sqrt(sig2))
x <- c(0, cumsum(x))
d <- as.data.frame(cbind(t,x))

plot(t, x, type = "l", ylim = c(-2, 2))

d |> ggplot(aes(t,x)) +
  geom_point(color = 'green') +
  geom_line(color = 'red')

5 Mô phỏng doanh thu của một quán cà phê

Trong ví dụ này chúng ta sẽ mô phỏng về doanh thu cho một quá cà phê nhỏ với những thông tin cơ bản sau:

  1. Diện tích khoảng \(40m^2\) bao gồm cả khu pha chế (khoảng \(10m^2\)).
  2. Đặt được 12 bàn, mỗi bàn 2 ghế.
  3. Quán bán những món sau: a. Cà phê đen - 20k/ly. b. Cà phê sữa đá - 25k/ly. c. Sữa chua - 15k/hủ. d. Sữa chua hạt đát - 20k/hủ.
  4. Giả sử về phân phối
    1. Số ly cà phê đen bán được trong ngày có phân phối đều trong khoảng 15 - 30.
    2. Số ly cà phê sữa đá có phân phối tam giác a = 20, b = 40, c = 35.
    3. Số hủ sữa chua bán được trong ngày có phân phối Poisson với tham số \(\lambda = 6\).
    4. Số ly sữa chua hạt đát bán được trong ngày có phân phối chuẩn với kỳ vọng \(\mu = 10\) và phương sai \(\sigma^2 = 5\).
library(EnvStats)
## 
## Attaching package: 'EnvStats'
## The following objects are masked from 'package:stats':
## 
##     predict, predict.lm
## The following object is masked from 'package:base':
## 
##     print.default
n <- 100

cfd <- runif(n, min = 15, max = 30)
cfd <- round(cfd,0)

cfsd <- rtri(n, min = 15, max = 35, mode = 25)
cfsd <- round(cfsd, digits = 0)

sc <- rpois(n, lambda = 6)

schd <- rnorm(n, mean = 10, sd = sqrt(5))
schd <- round(schd, digits = 0)

doanhthu <- cfd*20000 + cfsd*25000 + sc*15000 + schd*20000

hist(doanhthu)

doanhthuthap <- doanhthu[doanhthu < 1500000]
doanhthuvua <- doanhthu[doanhthu >= 1500000 & doanhthu <= 2200000]

doanhthucao <- doanhthu[doanhthu > 2200000]

table(cut(doanhthu,breaks = 3))
## 
## (1.03e+06,1.25e+06] (1.25e+06,1.47e+06] (1.47e+06,1.69e+06] 
##                  20                  50                  30
table(cut(doanhthu,3, labels = c('thap','vua','cao')))
## 
## thap  vua  cao 
##   20   50   30
table(cut(doanhthu/1000,3))
## 
## (1.03e+03,1.25e+03] (1.25e+03,1.47e+03] (1.47e+03,1.69e+03] 
##                  20                  50                  30
length(doanhthuthap)/length(doanhthu)
## [1] 0.78
length(doanhthuvua)/length(doanhthu)
## [1] 0.22
length(doanhthucao)/length(doanhthu)
## [1] 0