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)
- 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.
- 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.
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')
library(TeachingDemos)
twodice<- dice(200, ndice = 2, plot.it = T)
twodice$sum <- twodice$Red + twodice$Green
ggplot(twodice, aes(sum)) +
geom_bar(fill = 'red')
tmp <- runif(1000, 5, 10)
hist(tmp)
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')
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')
tmp <- rnorm(100)
hist(tmp)
tmp <- as.data.frame(tmp)
ggplot(data = tmp, aes(tmp)) +
geom_density(color = 'green')
tmp <- rexp(1000,0.5)
hist(tmp)
tmp <- as.data.frame(tmp)
ggplot(data = tmp, aes(tmp)) +
geom_density(color = 'green')
Để 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:
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.
Đâ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
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).
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')
Đị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:
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')
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:
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')
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:
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