1 Bài tập tuần 5

Chạy mô hình với các phân phối khác nhau của biến đầu vào.

1.1 Dữ liệu đầu vào

Dữ liệu gồm 52 quan sát từ quý 1 năm 2010 đến quý 4 năm 2022 gồm 6 biến:

CN: Công nghiệp
KK: Khai khoáng
CNCBCT: Công nghiệp chế biến, chế tạo
SXPP: Sản xuất và phân phối điện, khí đốt, nước nóng, hơi nước và điều hòa không khí
CC: Cung cấp nước, hoạt động quản lý và xử lý rác thải, nước thải
XD: Xây dựng

setwd("D:/HỌC TẬP/Mô phỏng ngẫu nhiên")
library(xlsx)
library(readxl)

## Warning: package 'readxl' was built under R version 4.2.3

library(data.table)

## Warning: package 'data.table' was built under R version 4.2.3

d <- read.xlsx("GDP_MPNN.xlsx",1)
data.table(d)

##         NA.       CN       KK   CNCBCT     SXPP        CC       XD
##  1: Q1 2010 27.80233 7.276075 17.72502 2.396061 0.4051811 5.518612
##  2: Q2 2010 26.75083 6.744417 17.06550 2.518298 0.4226227 6.085374
##  3: Q3 2010 26.97934 6.655979 17.36618 2.516098 0.4410775 6.188818
##  4: Q4 2010 25.52245 6.524900 16.39179 2.153805 0.4519643 7.249104
##  5: Q1 2011 29.69929 7.859633 19.33079 2.136214 0.3726605 5.086161
##  6: Q2 2011 29.79878 8.670915 18.52638 2.218944 0.3825372 5.293760
##  7: Q3 2011 29.26935 6.657590 19.74397 2.433711 0.4340879 5.254924
##  8: Q4 2011 27.80478 8.007136 17.24470 2.072588 0.4803541 6.158234
##  9: Q1 2012 31.43137 7.133887 21.70656 2.191048 0.3998756 4.703270
## 10: Q2 2012 31.89412 8.970369 20.28061 2.242091 0.4010558 5.225348
## 11: Q3 2012 29.68559 6.450603 20.17093 2.618891 0.4451643 5.156563
## 12: Q4 2012 29.25116 7.527098 19.01356 2.198952 0.5115411 6.119646
## 13: Q1 2013 30.76942 6.487755 21.44801 2.416436 0.4172150 4.609639
## 14: Q2 2013 31.39315 7.787355 20.63107 2.552621 0.4221004 5.043063
## 15: Q3 2013 30.67386 6.215322 21.35042 2.661172 0.4469455 4.860839
## 16: Q4 2013 29.15255 6.872262 19.37957 2.400505 0.5002216 5.823399
## 17: Q1 2014 30.51610 6.610706 20.73667 2.739031 0.4296976 4.609448
## 18: Q2 2014 31.05277 7.409740 20.28529 2.926901 0.4308440 4.933733
## 19: Q3 2014 30.47767 6.060404 21.04597 2.924489 0.4468084 4.947560
## 20: Q4 2014 29.00402 6.439169 19.43414 2.613998 0.5167128 5.667065
## 21: Q1 2015 29.23257 4.435080 21.35291 2.986814 0.4577655 4.749028
## 22: Q2 2015 29.45609 4.842444 20.90392 3.249980 0.4597477 5.185735
## 23: Q3 2015 29.64403 4.313458 21.52479 3.328762 0.4770228 5.255187
## 24: Q4 2015 27.01316 3.453826 20.07873 2.928111 0.5524936 6.524618
## 25: Q1 2016 29.07729 3.553725 21.68888 3.351740 0.4829513 4.734969
## 26: Q2 2016 29.03322 3.605723 21.55285 3.412853 0.4617858 5.306063
## 27: Q3 2016 29.37279 3.330322 22.16920 3.394086 0.4791873 5.318349
## 28: Q4 2016 26.97445 2.881695 20.57187 2.973055 0.5478269 6.631873
## 29: Q1 2017 29.17387 3.735911 21.48111 3.489877 0.4669709 4.711940
## 30: Q2 2017 29.56666 3.421134 22.23031 3.446898 0.4683159 5.381099
## 31: Q3 2017 30.68682 3.192556 23.42899 3.588540 0.4767325 5.427479
## 32: Q4 2017 29.72389 2.842052 23.19806 3.151522 0.5322520 6.659747
## 33: Q1 2018 30.47793 4.029690 22.41827 3.557580 0.4723820 4.800841
## 34: Q2 2018 30.59450 3.601389 22.81979 3.690668 0.4826517 5.609320
## 35: Q3 2018 31.93339 3.395256 24.61641 3.427264 0.4944588 5.638618
## 36: Q4 2018 30.31750 3.105171 23.47818 3.197829 0.5363122 6.626978
## 37: Q1 2019 30.68016 3.463981 22.83401 3.900762 0.4814111 4.839561
## 38: Q2 2019 30.96742 3.146650 23.46250 3.871339 0.4869268 5.669476
## 39: Q3 2019 32.18733 2.937786 24.93589 3.806626 0.5070302 5.760281
## 40: Q4 2019 29.97142 2.528621 23.79019 3.120956 0.5316565 6.978980
## 41: Q1 2020 30.34627 2.901958 22.88656 4.074066 0.4836864 4.908820
## 42: Q2 2020 30.11332 2.624843 22.77702 4.204425 0.5070250 5.918180
## 43: Q3 2020 31.76856 2.190370 24.95579 4.102598 0.5197984 5.890936
## 44: Q4 2020 30.67356 1.974749 24.88840 3.310571 0.4998360 7.065943
## 45: Q1 2021 30.57364 2.821993 23.22615 4.048846 0.4766484 4.908126
## 46: Q2 2021 31.35697 2.641854 23.97584 4.217828 0.5214420 5.870273
## 47: Q3 2021 32.86119 2.458709 25.44723 4.423421 0.5318323 6.055922
## 48: Q4 2021 31.28222 1.920831 25.62377 3.263171 0.4744490 6.854210
## 49: Q1 2022 32.02314 3.092619 24.21416 4.242842 0.4735109 5.054937
## 50: Q2 2022 32.81631 3.197458 24.98466 4.130962 0.5032227 6.050760
## 51: Q3 2022 32.78072 2.844548 25.02685 4.407470 0.5018547 6.490710
## 52: Q4 2022 30.79383 2.243662 24.76535 3.311026 0.4737833 6.989623
##         NA.       CN       KK   CNCBCT     SXPP        CC       XD

1.2 Mô phỏng

library(fBasics)

## Warning: package 'fBasics' was built under R version 4.2.3

basicStats(d$CN)

##                 X..d.CN
## nobs          52.000000
## NAs            0.000000
## Minimum       25.522454
## Maximum       32.861191
## 1. Quartile   29.246508
## 3. Quartile   30.988757
## Mean          30.046214
## Median        30.331883
## Sum         1562.403123
## SE Mean        0.223032
## LCL Mean      29.598458
## UCL Mean      30.493970
## Variance       2.586659
## Stdev          1.608309
## Skewness      -0.602237
## Kurtosis       0.213309

library(fBasics)
basicStats(d$CNCBCT)

##             X..d.CNCBCT
## nobs          52.000000
## NAs            0.000000
## Minimum       16.391785
## Maximum       25.623769
## 1. Quartile   20.253186
## 3. Quartile   23.466423
## Mean          21.734342
## Median        21.620866
## Sum         1130.185762
## SE Mean        0.331439
## LCL Mean      21.068951
## UCL Mean      22.399732
## Variance       5.712282
## Stdev          2.390038
## Skewness      -0.309763
## Kurtosis      -0.730576

library(fBasics)
basicStats(d$CC)

##               X..d.CC
## nobs        52.000000
## NAs          0.000000
## Minimum      0.372661
## Maximum      0.552494
## 1. Quartile  0.446397
## 3. Quartile  0.502197
## Mean         0.472724
## Median       0.476690
## Sum         24.581638
## SE Mean      0.005888
## LCL Mean     0.460903
## UCL Mean     0.484544
## Variance     0.001803
## Stdev        0.042458
## Skewness    -0.294518
## Kurtosis    -0.494198

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.2.3

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:data.table':
## 
##     between, first, last

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

n_bins <- 10
breaks <- quantile(d$KK, probs = seq(0, 1, length.out = n_bins + 1), type = 1, na.rm = TRUE)

kkg <- sapply(2:length(breaks), function(i) sample(d$KK[d$KK >= breaks[i - 1] & d$KK <= breaks[i]], 10, replace = TRUE))

K <- table(kkg)

Kk <- as.data.frame(K) %>%
  mutate(Tan_suat = Freq / sum(Freq))
Kk

##                 kkg Freq Tan_suat
## 1  1.97474913117254    1     0.01
## 2  2.19037016484278    2     0.02
## 3  2.24366243785466    2     0.02
## 4  2.45870893795037    2     0.02
## 5  2.52862056159668    5     0.05
## 6  2.62484344167525    1     0.01
## 7  2.64185419294828    1     0.01
## 8  2.82199267200874    1     0.01
## 9  2.84205178770792    4     0.04
## 10 2.84454830489354    4     0.04
## 11 2.88169508218852    2     0.02
## 12 2.90195809364746    2     0.02
## 13 2.93778579138361    2     0.02
## 14 3.09261934876639    1     0.01
## 15 3.10517103176282    1     0.01
## 16 3.14664992520624    4     0.04
## 17 3.19255628017669    1     0.01
## 18 3.33032150313402    2     0.02
## 19 3.39525603420142    6     0.06
## 20 3.46398112064342    1     0.01
## 21 3.55372534567234    4     0.04
## 22 3.60138947622404    2     0.02
## 23 3.60572294915913    1     0.01
## 24 3.73591123782747    2     0.02
## 25 4.02968977795517    2     0.02
## 26 4.31345833990272    3     0.03
## 27 4.84244445505115    1     0.01
## 28 6.06040444761687    3     0.03
## 29 6.21532187544485    2     0.02
## 30 6.43916904657199    2     0.02
## 31 6.45060306543078    1     0.01
## 32 6.48775540590923    2     0.02
## 33 6.52489964192492    1     0.01
## 34 6.61070647937773    3     0.03
## 35 6.65597943580786    4     0.04
## 36 6.65758990825106    1     0.01
## 37 6.74441653045645    3     0.03
## 38 6.87226183056873    1     0.01
## 39 7.13388654804438    2     0.02
## 40 7.27607503110032    2     0.02
## 41 7.40973976407404    1     0.01
## 42 7.52709766389301    4     0.04
## 43 7.85963287752904    3     0.03
## 44 8.67091503086552    2     0.02
## 45 8.97036938810019    3     0.03

library(dplyr)
n_bins <- 10
breaks <- quantile(d$SXPP, probs = seq(0, 1, length.out = n_bins + 1), type = 1, na.rm = TRUE)

sx <- sapply(2:length(breaks), function(i) sample(d$SXPP[d$SXPP >= breaks[i - 1] & d$SXPP <= breaks[i]], 10, replace = TRUE))

S <- table(sx)

Sx <- as.data.frame(S) %>%
  mutate(Tan_suat = Freq / sum(Freq))

Sx

##                  sx Freq Tan_suat
## 1  2.15380480044362    1     0.01
## 2  2.19104826362534    1     0.01
## 3  2.19895231629668    6     0.06
## 4  2.21894399118426    4     0.04
## 5  2.24209081149811    1     0.01
## 6  2.39606075386324    2     0.02
## 7  2.40050481276904    4     0.04
## 8  2.43371067088471    3     0.03
## 9  2.51609754688072    4     0.04
## 10  2.5182978121541    1     0.01
## 11 2.55262124062935    2     0.02
## 12 2.61889137813463    6     0.06
## 13 2.66117191348481    2     0.02
## 14 2.92811096947551    4     0.04
## 15 2.97305489385597    2     0.02
## 16 2.98681368388827    1     0.01
## 17 3.12095626761687    3     0.03
## 18 3.15152237782226    1     0.01
## 19 3.19782920657643    4     0.04
## 20 3.26317077436855    2     0.02
## 21 3.31057120569083    2     0.02
## 22 3.32876151068537    2     0.02
## 23 3.35173993702559    4     0.04
## 24 3.39408641098406    2     0.02
## 25 3.41285302107156    1     0.01
## 26 3.42726380961728    1     0.01
## 27 3.44689779613138    3     0.03
## 28 3.48987650424477    1     0.01
## 29 3.55758040201798    2     0.02
## 30 3.58853953069416    1     0.01
## 31 3.69066752906663    4     0.04
## 32 3.80662598727375    3     0.03
## 33 3.90076233030699    3     0.03
## 34 4.07406596990505    1     0.01
## 35 4.10259846919176    4     0.04
## 36 4.13096207790727    2     0.02
## 37 4.20442527948733    5     0.05
## 38 4.21782829225471    2     0.02
## 39 4.40747024247083    1     0.01
## 40 4.42342102715324    2     0.02

library(dplyr)
n_bins <- 10
breaks <- quantile(d$XD, probs = seq(0, 1, length.out = n_bins + 1), type = 1, na.rm = TRUE)

xaydung <- sapply(2:length(breaks), function(i) sample(d$XD[d$XD >= breaks[i - 1] & d$XD <= breaks[i]], 10, replace = TRUE))

X <- table(xaydung)

Xd <- as.data.frame(X) %>%
  mutate(Tan_suat = Freq / sum(Freq))
Xd

##             xaydung Freq Tan_suat
## 1  4.60944800824035    3     0.03
## 2  4.60963901748547    1     0.01
## 3  4.70326981177233    1     0.01
## 4  4.71193988910874    1     0.01
## 5  4.73496909417757    2     0.02
## 6  4.74902765309405    3     0.03
## 7  4.83956113499323    2     0.02
## 8  4.86083870371692    2     0.02
## 9  4.90812640376935    2     0.02
## 10 4.90881998434628    4     0.04
## 11 4.93373262353868    1     0.01
## 12 4.94755979332598    4     0.04
## 13 5.04306260223245    2     0.02
## 14  5.0549374862877    2     0.02
## 15 5.08616054360378    1     0.01
## 16 5.15656293376329    3     0.03
## 17 5.18573466785971    3     0.03
## 18 5.22534845930231    1     0.01
## 19 5.25492422210561    1     0.01
## 20 5.25518659722746    3     0.03
## 21 5.29376030400241    1     0.01
## 22 5.30606260262954    6     0.06
## 23 5.31834866586471    1     0.01
## 24  5.5186119455827    1     0.01
## 25 5.63861799100156    4     0.04
## 26 5.66947559458313    4     0.04
## 27 5.76028099735905    3     0.03
## 28 5.82339872880876    1     0.01
## 29 5.87027302648059    4     0.04
## 30 5.89093646773286    1     0.01
## 31 5.91818006317723    1     0.01
## 32 6.05075983543772    2     0.02
## 33 6.05592220052127    4     0.04
## 34 6.08537412367017    1     0.01
## 35 6.11964647587315    2     0.02
## 36 6.15823362860776    1     0.01
## 37 6.18881802299602    5     0.05
## 38 6.62697776316361    2     0.02
## 39 6.63187296744266    4     0.04
## 40 6.65974654307538    3     0.03
## 41   6.854209648257    1     0.01
## 42 6.97898022102724    1     0.01
## 43 6.98962291669853    1     0.01
## 44 7.06594290519492    3     0.03
## 45 7.24910423600482    1     0.01

# Công nghiệp 
cn <- rnorm(10000, mean = 30.046214, sd = 1.608309)

# Khai khoáng 
prob1 <- Kk$Tan_suat
khaikhoang1 <- sample(Kk$kkg, 10000, replace = TRUE , prob = prob1)
khaikhoang2 <-  as.character(khaikhoang1)
kk <- as.numeric(khaikhoang2)

# Công nghiệp chế biến, chế tạo 
cb <- rnorm(10000, mean = 21.734342, sd = 2.390038)

# Sản xuất và phân phối điện, khí đốt, nước nóng, hơi nước và điều hòa không khí 
prob2 <- Sx$Tan_suat
sanxuat1 <- sample(Sx$sx, 10000, replace = TRUE , prob = prob2)
sanxuat2 <-  as.character(sanxuat1)
sxpp <- as.numeric(sanxuat2)

# Cung cấp nước, hoạt động quản lý và xử lý rác thải, nước thải 
cc <- rnorm(10000, mean = 0.472724, sd = 0.042458)


# Xây dựng
prob3 <- Xd$Tan_suat
xaydung1 <- sample(Xd$xaydung, 10000, replace = TRUE , prob = prob3)
xaydung2 <-  as.character(xaydung1)
xd <- as.numeric(xaydung2)


GDP <- cn + kk + cb + sxpp + cc + xd 

hist(GDP)

summary(GDP)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   52.45   63.17   65.73   65.74   68.28   80.37

GDPthap <- GDP[GDP < 60]
GDPvua <- GDP[GDP >= 60 & GDP <= 70]

GDPcao <- GDP[GDP > 70]

table(cut(GDP,breaks = 3))

## 
## (52.4,61.8] (61.8,71.1] (71.1,80.4] 
##        1433        7793         774

table(cut(GDP,3, labels = c('thap','vua','cao')))

## 
## thap  vua  cao 
## 1433 7793  774

length(GDPthap)/length(GDP)

## [1] 0.0557

length(GDPvua)/length(GDP)

## [1] 0.8216

length(GDPcao)/length(GDP)

## [1] 0.1227

2 Bài tập tuần 4

Xây dụng mô hình cho đối tượng cần mô phỏng và giải thích mô hình.

2.1 Mô hình mô phỏng

Tổng sản phẩm trong nước là giá trị sản phẩm vật chất và dịch vụ cuối cùng được tạo ra của nền kinh tế trong một khoảng thời gian nhất định (quý, năm). Điều này có nghĩa trong GDP không tính các giá trị sản phẩm vật chất và dịch vụ đã sử dụng ở các khâu trung gian trong quá trình sản xuất tạo ra sản phẩm. GDP biểu thị kết quả sản xuất do các đơn vị thường trú tạo ra trong lãnh thổ kinh tế của một quốc gia.

Theo tiêu chuẩn phân loại ngành của Liên Hợp Quốc, có thể phân loại thành 3 nhóm ngành lớn:

Khu vực I: Nông nghiệp, lâm nghiệp và ngư nghiệp.
Khu vực II: Công nghiệp, xây dựng
Khu vực III: Dịch vụ
Cơ cấu ngành công nghiệp nước ta tương đối đa dạng bao gồm: Công nghiệp; khai khoáng; công nghiệp chế biến, chế tạo; sản xuất và phân phối điện, khí đốt, nước nóng, hơi nước và điều hòa không khí; Cung cấp nước, hoạt động quản lý và xử lý rác thải, nước thải; Xây dựng.

Tổng sản phẩm trong nước theo giá hiện hành phân theo nhóm ngành Công nghiệp và xây dựng, công thức tính:

  Tổng sản phẩm trong nước nhóm ngành CNXD  =   Công nghiệp + Khai khoáng + Công nghiệp chế biến, chế tạo + Sản xuất và phân phối điện, khí đốt, nước nóng, hơi nước và điều hòa không khí +  Cung cấp nước, hoạt động quản lý và xử lý rác thải, nước thải + Xây dựng.

3 Bài tập tuần 3:

Xác định phân phối cho các biến (ngẫu nhiên) đầu vào, giải thích.

3.1 Biến Công nghiệp

3.1.1 Đồ thị

CN <- d$CN
KK <- d$KK
CNCBCT <- d$CNCBCT
SXPP <- d$SXPP
CC <- d$CC
XD <- d$XD

hist(CN, col = 'pink')

3.1.2 Kiểm định phân phối

Kiểm định Shapiro-Wilk Test

Giả thuyết:

\(H_0\): CN tuân theo phân phối chuẩn

\(H_1\): CN không tuân theo phân phối chuẩn

C <- as.data.frame(CN)
shapiro.test(CN)

## 
##  Shapiro-Wilk normality test
## 
## data:  CN
## W = 0.9569, p-value = 0.05749

Kết quả kiểm định cho thấy p- value = 0.05749 > 0.05 nên chưa đủ cơ sở bác bỏ \(H_0\). Vậy CN tuân theo phân phối chuẩn.

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.2.3

d |> ggplot() +
  geom_qq(aes(sample = CN), color = 'red') +
  geom_qq_line(aes(sample = CN), color = 'green')

3.2 Biến Khai khoáng

3.2.1 Đồ thị

library(fBasics)
basicStats(KK)

##                     KK
## nobs         52.000000
## NAs           0.000000
## Minimum       1.920831
## Maximum       8.970369
## 1. Quartile   2.928829
## 3. Quartile   6.622025
## Mean          4.674834
## Median        3.603556
## Sum         243.091381
## SE Mean       0.288557
## LCL Mean      4.095532
## UCL Mean      5.254137
## Variance      4.329790
## Stdev         2.080815
## Skewness      0.457590
## Kurtosis     -1.327147

hist(KK, col = 'pink')

3.2.2 Kiểm định phân phối

Kiểm định Shapiro-Wilk Test

Giả thuyết:

\(H_0\): KK tuân theo phân phối chuẩn

\(H_1\): KK không tuân theo phân phối chuẩn

shapiro.test(KK)

## 
##  Shapiro-Wilk normality test
## 
## data:  KK
## W = 0.88306, p-value = 0.0001037

Kết quả kiểm định cho thấy p- value = 0.0001037 < 0.05 nên ta bác bỏ \(H_0\). Vậy KK không tuân theo phân phối chuẩn.

library(ggplot2)
d |> ggplot() +
  geom_qq(aes(sample = KK), color = 'red') +
  geom_qq_line(aes(sample = KK), color = 'green')

Kiểm định Kolmogorov-Smirnov (KS)

Giả thuyết

\(H_0\): KK tuân theo phân phối mũ

\(H_1\): KK không tuân theo phân phối mũ

ks.test(KK, y = 'pexp')

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  KK
## D = 0.85351, p-value < 2.2e-16
## alternative hypothesis: two-sided

Kết quả kiểm định cho thấy p- value < 0.05 nên ta bác bỏ \(H_0\). Vậy KK không tuân theo phân phối mũ.

ks.test(KK, y = 'plnorm')

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  KK
## D = 0.74504, p-value < 2.2e-16
## alternative hypothesis: two-sided

Kết quả kiểm định cho thấy p- value < 0.05 nên ta bác bỏ \(H_0\). Vậy KK không tuân theo phân phối loga chuẩn.

3.3 Biến Công nghiệp chế biến, chế tạo

3.3.1 Đồ thị

hist(CNCBCT, col = 'pink')

3.3.2 Kiểm định phân phối

Kiểm định Shapiro-Wilk Test

Giả thuyết:

\(H_0\): CNCBCT tuân theo phân phối chuẩn

\(H_1\): CNCBCT không tuân theo phân phối chuẩn

B <- as.data.frame(CNCBCT)
shapiro.test(CNCBCT)

## 
##  Shapiro-Wilk normality test
## 
## data:  CNCBCT
## W = 0.96832, p-value = 0.1791

Kết quả kiểm định cho thấy p- value = 0.1791 > 0.05 nên chưa đủ cơ sở bác bỏ \(H_0\). Vậy CNCBCT tuân theo phân phối chuẩn.

library(ggplot2)
d |> ggplot() +
  geom_qq(aes(sample = CNCBCT), color = 'red') +
  geom_qq_line(aes(sample = CNCBCT), color = 'green')

3.4 Biến SXPP

3.4.1 Đồ thị

library(fBasics)
basicStats(SXPP)

##                   SXPP
## nobs         52.000000
## NAs           0.000000
## Minimum       2.072588
## Maximum       4.423421
## 1. Quartile   2.544040
## 3. Quartile   3.614072
## Mean          3.164314
## Median        3.223905
## Sum         164.544342
## SE Mean       0.095360
## LCL Mean      2.972872
## UCL Mean      3.355757
## Variance      0.472860
## Stdev         0.687648
## Skewness      0.125870
## Kurtosis     -1.158070

hist(SXPP, col = 'pink')

3.4.2 Kiểm định phân phối

Kiểm định Shapiro-Wilk Test

Giả thuyết:

\(H_0\): SXPP tuân theo phân phối chuẩn

\(H_1\): SXPP không tuân theo phân phối chuẩn

shapiro.test(SXPP)

## 
##  Shapiro-Wilk normality test
## 
## data:  SXPP
## W = 0.9534, p-value = 0.04065

` Kết quả kiểm định cho thấy p- value = 0.04065 < 0.05 nên ta bác bỏ \(H_0\). Vậy SXPP không tuân theo phân phối chuẩn.

library(ggplot2)
d |> ggplot() +
  geom_qq(aes(sample = SXPP), color = 'red') +
  geom_qq_line(aes(sample = SXPP), color = 'green')

Kiểm định Kolmogorov-Smirnov (KS)

Giả thuyết

\(H_0\): SXPP tuân theo phân phối mũ

\(H_1\): SXPP không tuân theo phân phối mũ

ks.test(SXPP, y = 'pexp')

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  SXPP
## D = 0.87414, p-value < 2.2e-16
## alternative hypothesis: two-sided

Kết quả kiểm định cho thấy p- value < 0.05 nên ta bác bỏ \(H_0\). Vậy SXPP không tuân theo phân phối mũ.

ks.test(SXPP, y = 'plnorm')

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  SXPP
## D = 0.76694, p-value < 2.2e-16
## alternative hypothesis: two-sided

Kết quả kiểm định cho thấy p- value < 0.05 nên ta bác bỏ \(H_0\). Vậy SXPP không tuân theo phân phối loga chuẩn.

3.5 Biến Cung cấp

3.5.1 Đồ thị

hist(CC, col = 'pink')

3.5.2 Kiểm định phân phối

Kiểm định Shapiro-Wilk Test

Giả thuyết:

\(H_0\): CC tuân theo phân phối chuẩn

\(H_1\): CC không tuân theo phân phối chuẩn

Cc <- as.data.frame(CC)
shapiro.test(CC)

## 
##  Shapiro-Wilk normality test
## 
## data:  CC
## W = 0.98178, p-value = 0.6039

Kết quả kiểm định cho thấy p- value = 0.6039 > 0.05 nên chưa đủ cơ sở bác bỏ \(H_0\). Vậy CC tuân theo phân phối chuẩn.

library(ggplot2)
d |> ggplot() +
  geom_qq(aes(sample = CC), color = 'red') +
  geom_qq_line(aes(sample = CC), color = 'green')

3.6 Biến Xây dựng

3.6.1 Đồ thị

library(fBasics)
basicStats(XD)

##                     XD
## nobs         52.000000
## NAs           0.000000
## Minimum       4.609448
## Maximum       7.249104
## 1. Quartile   5.019187
## 3. Quartile   6.093942
## Mean          5.623137
## Median        5.473046
## Sum         292.403139
## SE Mean       0.102082
## LCL Mean      5.418199
## UCL Mean      5.828076
## Variance      0.541879
## Stdev         0.736124
## Skewness      0.515979
## Kurtosis     -0.863092

hist(XD, col = 'pink')

3.6.2 Kiểm định phân phối

Kiểm định Shapiro-Wilk Test

Giả thuyết:

\(H_0\): XD tuân theo phân phối chuẩn

\(H_1\): XD không tuân theo phân phối chuẩn

shapiro.test(XD)

## 
##  Shapiro-Wilk normality test
## 
## data:  XD
## W = 0.93825, p-value = 0.009511

Kết quả kiểm định cho thấy p- value = 0.002291 < 0.05 nên ta bác bỏ \(H_0\). Vậy XD không tuân theo phân phối chuẩn.

library(ggplot2)
d |> ggplot() +
  geom_qq(aes(sample = XD), color = 'red') +
  geom_qq_line(aes(sample = XD), color = 'green')

Kiểm định Kolmogorov-Smirnov (KS)

Giả thuyết

\(H_0\): XD tuân theo phân phối mũ

\(H_1\): XD không tuân theo phân phối mũ

ks.test(XD, y = 'pexp')

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  XD
## D = 0.99004, p-value < 2.2e-16
## alternative hypothesis: two-sided

Kết quả kiểm định cho thấy p- value < 0.05 nên ta bác bỏ \(H_0\). Vậy XD không tuân theo phân phối mũ.

ks.test(XD, y = 'plnorm')

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  XD
## D = 0.93676, p-value < 2.2e-16
## alternative hypothesis: two-sided

Kết quả kiểm định cho thấy p- value < 0.05 nên ta bác bỏ \(H_0\). Vậy XD không tuân theo phân phối loga chuẩn.

4 Bài tập tuần 2:

Chọn và giải thích về đối tượng (vấn đề) cần mô phỏng (có ít nhất 6 biến đầu vào)

4.1 Vấn đề cần mô phỏng

Đề tài: Mô phỏng các yếu tố cấu thành cơ cấu GDP nhóm ngành công nghiệp và xây dựng của Việt Nam.

Theo Bộ Công Thương, trong 10 năm (2011-2020), công nghiệp là ngành có tốc độ tăng trưởng cao nhất trong các ngành kinh tế quốc dân với đóng góp xấp xỉ 30% vào GDP và trở thành ngành xuất khẩu chủ lực của đất nước, góp phần đưa Việt Nam từ vị trí thứ 50 (năm 2010) lên vị trí thứ 22 (năm 2019) trong các quốc gia xuất khẩu lớn nhất thế giới.

Bộ Công Thương đặt ra mục tiêu đến năm 2030 Việt Nam là nước có công nghiệp hiện đại, thuộc nhóm quốc gia có năng lực cạnh tranh công nghiệp cao, nhóm 15 quốc gia xuất khẩu lớn nhất thế giới.

Nhận thấy tiềm năng tăng trưởng kinh tế Việt Nam thuộc thuộc nhiều vào nhóm ngành Công nghiệp và xây dựng, tôi sẽ đi mô phỏng các yếu tố cấu thành nhóm ngành này như: Công nghiệp; Khai khoáng; Công nghiệp chế biến, chế tạo; Sản xuất và phân phối ; Cung cấp ; xây dựng.

4.2 Dữ liệu đầu vào

Hệ thống phổ biến dữ liệu chung tăng cường (e-GDDS) Trang dữ liệu tóm tắt quốc gia (NSDP) – VIỆT NAM http://nsdp.gso.gov.vn/

Dữ liệu gồm 52 quan sát từ quý 1 năm 2010 đến quý 4 năm 2022 gồm 6 biến:

CN: Công nghiệp
KK: Khai khoáng
CNCBCT: Công nghiệp chế biến, chế tạo
SXPP: Sản xuất và phân phối điện, khí đốt, nước nóng, hơi nước và điều hòa không khí
CC: Cung cấp nước, hoạt động quản lý và xử lý rác thải, nước thải
XD: Xây dựng

5 BÀI VỀ NHÀ TUẦN 1:

Mô phỏng ít nhất 5 biến ngẫu nhiên (có phân phối xác suất khác nhau), mô phỏng, vẽ đồ thị, tính toàn các đặc trưng đo lường và giải thích ý nghĩa.

5.1 Phân phối đều

X ∼ Unif{1, 2, · · · , n} là mô hình đều rời rạc trên tập hợp số tự nhiên từ 1, 2, · · · , n. Hàm P(X = x) = 1/n được gọi là hàm xác suất khối lượng của biến X tại điểm x.

cpp1 <- runif(n = 1000,min = 8,max = 10) #Sinh 1000 biến cpp1 có phân phối đều trong khoảng [8,10]
hist(cpp1, main = "Uniform Distribution", xlab = "ccp1")

summary(cpp1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   8.009   8.532   9.050   9.026   9.505  10.000

5.2 Phân phối chuẩn chính tắc

Phân phối chuẩn tắc (standard normal distribution) là phân phối chuẩn với giá trị trung bình (μ) bằng 0 và độ lệch chuẩn (σ) bằng 1. Phân phối chuẩn còn được gọi là đường cong chuông (bell curve) vì đồ thị của mật độ xác suất có dạng chuông.

cpp2 <- rnorm(100) # Sinh 100 biến cpp2 có phân phối chuẩn trung bình 0 và độ lệch tiêu chuẩn 1
hist(cpp2, main = "Normal Distribution", xlab = "cpp2")

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

summary(cpp2)

##       cpp2        
##  Min.   :-1.6167  
##  1st Qu.:-0.4420  
##  Median : 0.2766  
##  Mean   : 0.2398  
##  3rd Qu.: 0.8226  
##  Max.   : 3.7247

5.3 Phân phối Poisson

Biến ngẫu nhiên X có phân phối Poisson là biến ngẫu nhiên dùng để mô tả cho số lần xảy ra của một sự việc/biến cố (event) mà chúng ta quan tâm xảy ra trong một khoảng thời gian hoặc không gian cho trước.

cpp3 <- rpois(n = 10,lambda = 1.2)
cpp3

##  [1] 1 3 0 3 2 2 2 1 0 0

cpp3 <- table(cpp3)
cpp3 <- as.data.frame(cpp3)
ggplot(cpp3,aes(cpp3,Freq)) + geom_col(fill='red') +
  geom_point() + geom_line(aes(as.integer(cpp3),Freq),color='green')

summary(cpp3)

##  cpp3       Freq    
##  0:1   Min.   :2.0  
##  1:1   1st Qu.:2.0  
##  2:1   Median :2.5  
##  3:1   Mean   :2.5  
##        3rd Qu.:3.0  
##        Max.   :3.0

5.4 Phân phối nhị thức

Nếu phép thử Bernoulli được thực hiện lặp lại n lần (độc lập nhau). Khi đó số phép thử thành công X trong n phép thử được thực hiện sẽ tuân theo mô hình nhị thức (binomial model) với hai tham số n và p, ký hiệu X ∼ Binomial(n, p), n ≥ 1, p ∈ (0, 1).

cpp4 <- rbinom(n = 100, size = 15, prob = 0.2)
cpp4

##   [1] 3 3 4 4 5 1 1 2 3 0 4 5 3 3 6 3 2 1 2 0 3 4 5 4 3 6 7 4 4 2 7 6 5 1 3 4 2
##  [38] 2 2 2 4 7 4 5 4 6 4 3 5 3 2 4 4 0 1 3 5 3 4 3 2 5 4 5 2 5 2 2 1 3 4 4 3 1
##  [75] 2 1 1 3 2 3 2 2 5 2 2 4 2 4 1 3 1 1 5 2 2 0 2 3 3 3

cpp4 <- table(cpp4)
cpp4 <- as.data.frame(cpp4)
ggplot(cpp4,aes(cpp4,Freq)) + geom_col(fill='red')

summary(cpp4)

##       cpp4        Freq     
##  0      :1   Min.   : 3.0  
##  1      :1   1st Qu.: 4.0  
##  2      :1   Median :12.0  
##  3      :1   Mean   :12.5  
##  4      :1   3rd Qu.:20.5  
##  5      :1   Max.   :23.0  
##  (Other):2

5.5 Phân phối siêu bội

Giả sử tổng thể có N phần tử, trong đó có n (n ≤ N) phần tử có dấu hiệu A và (N-n) phần tử có dấu hiệu không A. Lấy ngẫu nhiên một mẫu có kích thước k, k ≤ N. Số phần tử có dấu hiệu A trong mẫu được lấy ra, là một biến ngẫu nhiên có mô hình siêu bội, ký hiệu X ∼ Hypergeometric(N, n, k, x).

cpp5 <- rhyper(nn = 100,m = 100,n = 75,k = 25)
hist(cpp5, main = " Hypergeometric Distribution", xlab = "cpp5")

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

summary(cpp5)

##       cpp5      
##  Min.   : 7.00  
##  1st Qu.:13.00  
##  Median :14.00  
##  Mean   :14.24  
##  3rd Qu.:16.00  
##  Max.   :19.00

BTVN

Cao Phạm Hằng

2023-06-05

1 Bài tập tuần 5

1.1 Dữ liệu đầu vào

1.2 Mô phỏng

2 Bài tập tuần 4

2.1 Mô hình mô phỏng

3 Bài tập tuần 3:

3.1 Biến Công nghiệp

3.1.1 Đồ thị

3.1.2 Kiểm định phân phối

3.2 Biến Khai khoáng

3.2.1 Đồ thị

3.2.2 Kiểm định phân phối

3.3 Biến Công nghiệp chế biến, chế tạo

3.3.1 Đồ thị

3.3.2 Kiểm định phân phối

3.4 Biến SXPP

3.4.1 Đồ thị

3.4.2 Kiểm định phân phối

3.5 Biến Cung cấp

3.5.1 Đồ thị

3.5.2 Kiểm định phân phối

3.6 Biến Xây dựng

3.6.1 Đồ thị

3.6.2 Kiểm định phân phối

4 Bài tập tuần 2:

4.1 Vấn đề cần mô phỏng

4.2 Dữ liệu đầu vào

5 BÀI VỀ NHÀ TUẦN 1:

5.1 Phân phối đều

5.2 Phân phối chuẩn chính tắc

5.3 Phân phối Poisson

5.4 Phân phối nhị thức

5.5 Phân phối siêu bội