THNN2

1 Load các thư viện cần thiết

library(DT)           # Để hiển thị bảng dữ liệu
library(triangle)     # Để tạo phân phối tam giác
library(lubridate)    # Để xử lý ngày tháng
library(dplyr)        # Để sử dụng hàm case_when
library(knitr)        # Để tạo bảng
library(ggplot2)      # Để vẽ biểu đồ
library(patchwork)    # Để nối 2 biểu đồ
library(ggpubr)       # Cho các biểu đồ thống kê nâng cao
library(tidyverse)    # Tập hợp các package quan trọng cho data science
library(skimr)        # Xem tổng quan dữ liệu với thống kê mô tả
library(gt)           # Tạo bảng biểu đẹp và chuyên nghiệp
library(RColorBrewer) # Cung cấp các palette màu để vẽ biểu đồ đẹp
library(openxlsx)     # Xuất dữ liệu ra file exel

2 Biểu đồ các phân phối

2.1 Phân phối đều

# Tạo dữ liệu từ phân phối đều trong khoảng [5, 10]
ppDeu <- runif(1000, min = 5, max = 10)

# Vẽ biểu đồ histogram với trục Y là phần trăm
hist(ppDeu,
     breaks = 20,
     col = "lightgreen",
     freq = FALSE,  # biểu đồ theo mật độ
     main = "Phân phối Đều (Uniform)",
     xlab = "Giá trị",
     ylab = "Tần suất (%)")

# Chuyển trục Y từ mật độ sang phần trăm bằng cách nhân với 100
# (vì mật độ có tổng là 1, nên nhân 100 để thành phần trăm)
lines(density(ppDeu), col = "darkgreen", lwd = 2)

2.2 Phân phối chuẩn

# Tạo dữ liệu từ phân phối chuẩn (mean = 0, sd = 1)
ppChuan <- rnorm(1000, mean = 0, sd = 1)

# Vẽ biểu đồ histogram với trục Y là phần trăm
hist(ppChuan,
     breaks = 20,
     col = "skyblue",
     freq = FALSE,  # biểu đồ theo mật độ
     main = "Phân phối Chuẩn (Normal)",
     xlab = "Giá trị",
     ylab = "Tần suất (%)")

# Vẽ đường mật độ và nhân với 100 để chuyển sang phần trăm trực quan
lines(density(ppChuan), col = "blue", lwd = 2)

2.3 Phân phối Tam giác

# Tạo dữ liệu từ phân phối Tam giác (giá trị min = 0, mode = 5, max = 10)
ppTamGiac <- rtriangle(1000, a = 0, b = 10, c = 5)

# Vẽ biểu đồ histogram với trục Y là phần trăm
hist(ppTamGiac,
     breaks = 20,
     col = "lightblue",
     freq = FALSE,  # để biểu đồ chuẩn hóa
     main = "Phân phối Tam giác",
     xlab = "Giá trị",
     ylab = "Tần suất (%)")

# Chuyển tần suất sang phần trăm
lines(density(ppTamGiac), col = "blue", lwd = 2)

2.4 Phân phối rời rạc

ppRoiRac <- sample(c(1, 2, 3), size = 1000, replace = TRUE, prob = c(0.2, 0.5, 0.3))

table_pct <- prop.table(table(ppRoiRac)) * 100

barplot(table_pct,
        col = "salmon",
        main = "Phân phối Rời rạc",
        ylab = "Tần suất (%)",
        xlab = "Giá trị",
        ylim = c(0, max(table_pct) + 5))

3 Mô phỏng dữ liệu đơn hàng

# Khởi tạo seed để kết quả có thể tái lập
set.seed(123)

# 1. Tạo dữ liệu ngày trong năm 2025
start_date <- as.Date("2025-01-01")
end_date <- as.Date("2025-12-31")
dates <- seq.Date(start_date, end_date, by = "day")

# 2. Tạo mã đơn hàng theo tháng
generate_order_ids <- function(year, month) {
  days_in_month <- days_in_month(as.Date(paste(year, month, "01", sep = "-")))
  orders_per_day <- round(rtriangle(n = days_in_month, a = 10, b = 20, c = 15))
  
  order_ids <- character(0)
  for (day in 1:days_in_month) {
    day_orders <- sprintf("%s%02d%02d%03d", 
                          substr(year, 3, 4), 
                          month, 
                          day, 
                          1:orders_per_day[day])
    order_ids <- c(order_ids, day_orders)
  }
  return(order_ids)
}

# Tạo tất cả mã đơn hàng cho năm 2025
all_order_ids <- unlist(lapply(1:12, function(month) generate_order_ids(2025, month)))

# 3. Tạo loại sản phẩm theo phân phối rời rạc 
product_types <- sample(c("Decal", "Biểu mẫu", "Catalog"), 
                        size = length(all_order_ids), 
                        replace = TRUE, 
                        prob = c(0.5, 0.3, 0.2))

# Hàm làm tròn số lượng theo quy tắc thực tế
round_to_realistic <- function(x, product_type) {
  if (product_type == "Biểu mẫu") {
    round_to <- c(0)
  } else if (product_type == "Catalog") {
    round_to <- c(100, 500, 1000)
  } else { # Decal
    round_to <- c(100, 500, 1000)
  }
  
  differences <- sapply(round_to, function(r) abs(x - r))
  closest <- round_to[which.min(differences)]
  
  if (min(differences) > 0.2 * x) {
    return(round(x / 100) * 100)
  } else {
    return(closest)
  }
}

# 4. Tạo số lượng khách đặt cho từng loại sản phẩm
generate_product_data <- function(product_type) {
  if (product_type == "Biểu mẫu") {
    pages <- round(runif(1, 80, 100))  # Số trang mỗi cuốn biểu mẫu
    base_quantity <- rnorm(1, mean = 400, sd = 40)  # Số cuốn
    quantity <- round_to_realistic(base_quantity, product_type)
    return(list(quantity = quantity, pages_sheets = pages))
  } else if (product_type == "Catalog") {
    pages <- round(runif(1, 20, 40))  # Số trang mỗi cuốn catalog
    base_quantity <- rnorm(1, mean = 2000, sd = 200)  # Số cuốn
    quantity <- round_to_realistic(base_quantity, product_type)
    return(list(quantity = quantity, pages_sheets = pages))
  } else { # Decal
    sheets <- round(runif(1, 2000, 4000))  # Số tờ decal
    quantity <- round_to_realistic(sheets, product_type)
    return(list(quantity = quantity, pages_sheets = 0))  # Decal không có số trang
  }
}

# Tạo dữ liệu sản phẩm
product_data <- lapply(product_types, generate_product_data)
quantities <- sapply(product_data, function(x) x$quantity)
pages_sheets <- sapply(product_data, function(x) x$pages_sheets)

# 5. Tạo tỷ lệ lỗi 
error_rates <- numeric(length(product_types))

error_rates[product_types == "Biểu mẫu"] <-
  runif(sum(product_types == "Biểu mẫu"), 0.005, 0.01)

error_rates[product_types == "Catalog"] <-
  runif(sum(product_types == "Catalog"), 0.01, 0.02)

error_rates[product_types == "Decal"] <-
  runif(sum(product_types == "Decal"), 0.02, 0.03)

# 6. Tính số lượng lỗi
error_quantities <- round(quantities * error_rates)

# 7. Tạo đơn vị
units <- ifelse(product_types == "Decal", "Tờ", "Cuốn")

# 8. Tạo loại giấy
paper_types <- ifelse(product_types == "Decal", "CPWG-P0031",
                     ifelse(product_types == "Biểu mẫu", "Carbonless",
                            ifelse(runif(length(product_types)) > 0.5, "F100", "C100")))

# 9. Tạo ĐƠN GIÁ BÁN (giá bán cho khách hàng)
selling_prices <- case_when(
  paper_types == "CPWG-P0031" ~ 553,    # Giá bán mỗi tờ decal
  paper_types == "Carbonless" ~ 47,     # Giá bán mỗi trang biểu mẫu
  paper_types == "F100" ~ 30,           # Giá bán mỗi trang catalog
  paper_types == "C100" ~ 31            # Giá bán mỗi trang catalog
)

# 10. Tạo GIÁ VỐN (chi phí sản xuất)
cost_prices <- case_when(
  paper_types == "CPWG-P0031" ~ 388,    # Giá vốn mỗi tờ decal
  paper_types == "Carbonless" ~ 37,      # Giá vốn mỗi trang biểu mẫu
  paper_types == "F100" ~ 20,            # Giá vốn mỗi trang catalog
  paper_types == "C100" ~ 21             # Giá vốn mỗi trang catalog
)

# 11. Tính doanh thu theo từng loại sản phẩm (dùng giá bán)
calculate_revenue <- function(product_type, quantity, pages, selling_price) {
  if (product_type == "Decal") {
    # Decal: Số tờ × giá bán/tờ
    quantity * selling_price
  } else {
    # Biểu mẫu và Catalog: Số cuốn × số trang × giá bán/trang 
    quantity * pages * selling_price 
  }
}

revenue <- mapply(calculate_revenue, 
                 product_types, quantities, pages_sheets, selling_prices)

# 12. Tính chi phí biến động (dùng giá vốn)
calculate_variable_cost <- function(product_type, error_quantity, pages, cost_price) {
  if (product_type == "Decal") {
    # Decal: Số tờ lỗi × giá vốn/tờ
    error_quantity * cost_price
  } else {
    # Biểu mẫu và Catalog: Số cuốn lỗi × số trang × giá vốn/trang
    error_quantity * pages * cost_price
  }
}

variable_costs <- mapply(calculate_variable_cost, 
                        product_types, error_quantities, pages_sheets, cost_prices)

# Tạo data frame tổng hợp
daily_orders <- round(rtriangle(n = length(dates), a = 10, b = 20, c = 15))
data <- data.frame(
  Ngay = rep(dates, times = daily_orders)[1:length(revenue)],
  MaDonHang = all_order_ids[1:length(revenue)],
  LoaiSanPham = product_types[1:length(revenue)],
  SoLuongDatHang = quantities[1:length(revenue)],
  DonVi = units[1:length(revenue)],
  SoTrang = ifelse(product_types[1:length(revenue)] %in% c("Biểu mẫu", "Catalog"), 
                  pages_sheets[1:length(revenue)], 0),
  TyLeLoi = error_rates[1:length(revenue)],
  SoLuongLoi = error_quantities[1:length(revenue)],
  LoaiGiay = paper_types[1:length(revenue)],
  DonGiaBan = selling_prices[1:length(revenue)],
  GiaVon = cost_prices[1:length(revenue)],
  DoanhThu = revenue[1:length(revenue)],
  ChiPhiBienDong = variable_costs[1:length(revenue)]
)

# Sắp xếp lại theo ngày và làm sạch dữ liệu
data <- data[order(data$Ngay), ]
rownames(data) <- NULL

# Tạo bảng tổng hợp doanh thu và chi phí theo ngày
profit_data <- data %>% 
  group_by(Ngay) %>% 
  summarise(
    TongDoanhThu = sum(DoanhThu, na.rm = TRUE),
    TongChiPhiBienDong = sum(ChiPhiBienDong, na.rm = TRUE),
    ChiPhiCoDinh = 15000000,  # 15 triệu đồng/ngày
    .groups = "drop"
  ) %>% 
  mutate(
    LoiNhuan = TongDoanhThu - ChiPhiCoDinh - TongChiPhiBienDong
  )

# Tạo workbook mới
wb <- createWorkbook()

# Thêm sheet đầu tiên: dữ liệu gốc
addWorksheet(wb, "Data")
writeData(wb, sheet = "Data", data)

# Thêm sheet thứ hai: dữ liệu lợi nhuận
addWorksheet(wb, "Profit_Data")
writeData(wb, sheet = "Profit_Data", profit_data)

# Ghi file Excel ra đĩa
saveWorkbook(wb, "VuongNhuThuy2221000329_dulieu.xlsx", overwrite = TRUE)

datatable(data)

datatable(profit_data)

4 Thống kê mô tả các biến

4.1 Tổng quan

skim(data) %>% summary()

Data summary

Name	data
Number of rows	5483
Number of columns	13
_______________________
Column type frequency:
character	4
Date	1
numeric	8
________________________
Group variables	None

4.2 Thống kê mô tả biến Loại sản phẩm

data %>%
  count(LoaiSanPham) %>%
  mutate(Percent = round(n / sum(n) * 100, 2)) %>%
  kable(caption = "Bảng tần số và tần suất theo Loại sản phẩm")

Bảng tần số và tần suất theo Loại sản phẩm

LoaiSanPham	n	Percent
Biểu mẫu	1631	29.75
Catalog	1087	19.82
Decal	2765	50.43

# 1. Biểu đồ cột (KHÔNG hiển thị legend)
plot_bar <- ggplot(data, aes(x = LoaiSanPham, fill = LoaiSanPham)) +
  geom_bar() +
  geom_text(stat = 'count', aes(label = after_stat(count)),  position = position_stack(vjust = 0.5)) +
  labs(title = "SỐ LƯỢNG ĐƠN HÀNG",
       x = NULL,
       y = NULL) +
  theme_minimal() +
  scale_fill_brewer(palette = "Set4") +  # Dùng palette chung
  theme(legend.position = "none")  # Ẩn legend

# 2. Biểu đồ tròn (HIỂN THỊ legend)
plot_pie <- ggplot(data, aes(x = "", fill = LoaiSanPham)) +
  geom_bar(width = 0.5) +
  coord_polar("y") +
  geom_text(stat = "count",
            aes(label = scales::percent(after_stat(count)/sum(after_stat(count)))),
            position = position_stack(vjust = 0.5)) +
  labs(title = "TỶ LỆ PHÂN BỐ",
       x = NULL,
       y = NULL) +
  theme_void() +
  scale_fill_brewer(palette = "Set4")  # Dùng CÙNG palette với biểu đồ cột

# Kết hợp - chỉ hiển thị 1 legend
combined_plot <- plot_bar + plot_pie 

print(combined_plot)

4.3 Thống kê mô tả biến Loại giấy

data %>%
  count(LoaiGiay) %>%
  mutate(Percent = round(n / sum(n) * 100, 2)) %>%
  kable(caption = "Bảng tần số và tần suất theo Loại giấy")

Bảng tần số và tần suất theo Loại giấy

LoaiGiay	n	Percent
C100	522	9.52
CPWG-P0031	2765	50.43
Carbonless	1631	29.75
F100	565	10.30

# 1. Biểu đồ cột (KHÔNG hiển thị legend)
plot_bar <- ggplot(data, aes(x = LoaiGiay, fill = LoaiGiay)) +
  geom_bar() +
  geom_text(stat = 'count', aes(label = after_stat(count)),  position = position_stack(vjust = 0.5)) +
  labs(title = "SỐ LƯỢNG ĐƠN HÀNG",
       x = NULL,
       y = NULL) +
  theme_minimal() +
  scale_fill_brewer(palette = "Set4") +  # Dùng palette chung
  theme(legend.position = "none",  # Ẩn legend
  axis.text.x = element_text(size = 8, angle = 15, hjust = 1))
 
# 2. Biểu đồ tròn (HIỂN THỊ legend)
plot_pie <- ggplot(data, aes(x = "", fill = LoaiGiay)) +
  geom_bar(width = 0.5) +
  coord_polar("y") +
  geom_text(stat = "count",
            aes(label = scales::percent(after_stat(count)/sum(after_stat(count)))),
            position = position_stack(vjust = 0.5), size = 2.5) +
  labs(title = "TỶ LỆ PHÂN BỐ",
       x = NULL,
       y = NULL) +
  theme_void() +
  scale_fill_brewer(palette = "Set4")  # Dùng CÙNG palette với biểu đồ cột

# Kết hợp - chỉ hiển thị 1 legend
combined_plot <- plot_bar + plot_pie 

print(combined_plot)

4.4 Thống kê mô tả biến Loại sản phẩm và Số lượng đặt hàng

stats_table <- data %>%
  group_by(LoaiSanPham) %>%
  summarise(
    Mean = mean(SoLuongDatHang, na.rm = TRUE),
    SD = sd(SoLuongDatHang, na.rm = TRUE),
    Median = median(SoLuongDatHang, na.rm = TRUE),
    Min = min(SoLuongDatHang, na.rm = TRUE),
    Max = max(SoLuongDatHang, na.rm = TRUE)
  ) %>%
  mutate(across(where(is.numeric), ~ round(., 2)))  # Làm tròn 2 chữ số

# Hiển thị bảng
print(stats_table)

## # A tibble: 3 × 6
##   LoaiSanPham  Mean    SD Median   Min   Max
##   <chr>       <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1 Biểu mẫu     398.  45.3    400   300   500
## 2 Catalog     1994. 205.    2000  1300  2600
## 3 Decal       2989. 565.    3000  2000  4000

# 1. Boxplot - Phân bố số lượng đặt hàng theo loại sản phẩm
p1 <- ggplot(data, aes(x = LoaiSanPham, y = SoLuongDatHang, fill = LoaiSanPham)) +
  geom_boxplot(alpha = 0.7, outlier.shape = 19) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, color = "red") +
  labs(title = "PHÂN BỐ SỐ LƯỢNG ĐẶT HÀNG THEO LOẠI SP",
       x = NULL,
       y = "Số lượng đặt hàng") +
  theme_minimal() +
  scale_fill_brewer(palette = "Pastel1") +
  theme(legend.position = "none")

# 2. Violin plot + Jitter - Kết hợp phân bố và điểm dữ liệu
p2 <- ggplot(data, aes(x = LoaiSanPham, y = SoLuongDatHang, fill = LoaiSanPham)) +
  geom_violin(alpha = 0.5, trim = FALSE) +
  geom_jitter(width = 0.2, alpha = 0.3, size = 1.5) +
  labs(title = "MẬT ĐỘ PHÂN BỐ KÈM DATA POINTS",
       x = NULL,
       y = NULL) +
  theme_minimal() +
  scale_fill_brewer(palette = "Set2") +
  theme(legend.position = "none")

# 3. Barplot (Mean + SD)
p3 <- ggplot(data, aes(x = LoaiSanPham, y = SoLuongDatHang, fill = LoaiSanPham)) +
  stat_summary(fun = mean, geom = "bar", position = "dodge") +
  stat_summary(fun.data = mean_sdl, 
               geom = "errorbar", 
               width = 0.2,
               color = "gray30") +
  labs(title = "GIÁ TRỊ TRUNG BÌNH ± ĐỘ LỆCH CHUẨN",
       x = NULL,
       y = "Số lượng đặt hàng (Mean ± SD)") +
  theme_minimal() +
  scale_fill_brewer(palette = "Accent") +
  theme(legend.position = "none")

# 4. Phân bố tần suất (Histogram phân nhóm)
p4 <- ggplot(data, aes(x = SoLuongDatHang, fill = LoaiSanPham)) +
  geom_histogram(binwidth = 100, 
                 alpha = 0.7, 
                 position = "identity",
                 color = "white") +
  facet_wrap(~LoaiSanPham, ncol = 1) +
  labs(title = "PHÂN BỐ TẦN SỐ CỦA SỐ LƯỢNG ĐẶT HÀNG",
       x = "Số lượng đặt hàng",
       y = "Tần số") +
  theme_minimal() +
  scale_fill_brewer(palette = "Dark2")

# Kết hợp 4 biểu đồ
combined_plot <- (p1 + p2) / (p3 + p4) + 
  plot_layout(heights = c(1, 1.5)) +
  plot_annotation(title = "THỐNG KÊ MÔ TẢ: LOẠI SP vs SỐ LƯỢNG ĐẶT HÀNG",
                  theme = theme(plot.title = element_text(size = 16, hjust = 0.5, face = "bold")))

# Xuất file ảnh chất lượng cao (tùy chọn)
ggsave("LoaiSanPham_SoLuongDatHang.png", combined_plot, width = 12, height = 10, dpi = 300)

knitr::include_graphics("LoaiSanPham_SoLuongDatHang.png", error = FALSE)

data %>%
  filter(LoaiSanPham %in% c("Decal", "Biểu mẫu", "Catalog")) %>%
  count(LoaiSanPham, SoLuongDatHang) %>%
  group_by(LoaiSanPham) %>%
  gt() %>%
  tab_header(
    title = "BẢNG THỐNG KÊ TẦN SỐ CỦA SỐ LƯỢNG ĐẶT HÀNG",
    subtitle = "Phân theo Loại sản phẩm"
  ) %>%
  fmt_number(columns = n, decimals = 0) %>%
  cols_label(
    LoaiSanPham = "Loại sản phẩm",
    SoLuongDatHang = "Số lượng đặt hàng",
    n = "Tần số"
  )

BẢNG THỐNG KÊ TẦN SỐ CỦA SỐ LƯỢNG ĐẶT HÀNG
Phân theo Loại sản phẩm
Số lượng đặt hàng	Tần số
Biểu mẫu
300	183
400	1,296
500	152
Catalog
1300	2
1400	5
1500	9
1600	38
1700	74
1800	138
1900	177
2000	204
2100	182
2200	156
2300	65
2400	30
2500	5
2600	2
Decal
2000	53
2100	134
2200	140
2300	139
2400	151
2500	142
2600	133
2700	133
2800	161
2900	150
3000	150
3100	161
3200	132
3300	136
3400	140
3500	126
3600	133
3700	136
3800	135
3900	118
4000	62

# Lọc Decal và thống kê tần số
plot_freq <- data %>%
  filter(LoaiSanPham == "Decal") %>%
  ggplot(aes(x = SoLuongDatHang)) +
  geom_bar(fill = "steelblue") +
  geom_text(stat = 'count', aes(label = after_stat(count)), 
            vjust = -0.5, size = 3) +
  labs(title = "TẦN SỐ SỐ LƯỢNG ĐẶT HÀNG DECAL",
       x = "Số lượng đặt hàng",
       y = "Tần số") +
  theme_minimal()

print(plot_freq)

# Lọc Biểu mẫu và thống kê tần số
plot_freq <- data %>%
  filter(LoaiSanPham == "Biểu mẫu") %>%
  ggplot(aes(x = SoLuongDatHang)) +
  geom_bar(fill = "steelblue") +
  geom_text(stat = 'count', aes(label = after_stat(count)), 
            vjust = -0.5, size = 3) +
  labs(title = "TẦN SỐ SỐ LƯỢNG ĐẶT HÀNG BIỂU MẪU",
       x = "Số lượng đặt hàng",
       y = "Tần số") +
  theme_minimal()

print(plot_freq)

# Lọc Catalog và thống kê tần số
plot_freq <- data %>%
  filter(LoaiSanPham == "Catalog") %>%
  ggplot(aes(x = SoLuongDatHang)) +
  geom_bar(fill = "steelblue") +
  geom_text(stat = 'count', aes(label = after_stat(count)), 
            vjust = -0.5, size = 3) +
  labs(title = "TẦN SỐ SỐ LƯỢNG ĐẶT HÀNG CATALOG",
       x = "Số lượng đặt hàng",
       y = "Tần số") +
  theme_minimal()

print(plot_freq)

4.5 Thống kê mô tả biến Loại sản phẩm và Số trang

stats_table <- data %>%
  group_by(LoaiSanPham) %>%
  summarise(
    Mean = mean(SoTrang, na.rm = TRUE),
    SD = sd(SoTrang, na.rm = TRUE),
    Median = median(SoTrang, na.rm = TRUE),
    Min = min(SoTrang, na.rm = TRUE),
    Max = max(SoTrang, na.rm = TRUE)
  ) %>%
  mutate(across(where(is.numeric), ~ round(., 2)))  # Làm tròn 2 chữ số

# Hiển thị bảng
print(stats_table)

## # A tibble: 3 × 6
##   LoaiSanPham  Mean    SD Median   Min   Max
##   <chr>       <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1 Biểu mẫu     90.1  5.82     90    80   100
## 2 Catalog      30.2  5.75     30    20    40
## 3 Decal         0    0         0     0     0

# 1. Boxplot - Giữ nguyên
p1 <- ggplot(data, aes(x = LoaiSanPham, y = SoTrang, fill = LoaiSanPham)) +
  geom_boxplot(alpha = 0.7, outlier.shape = 19) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, color = "red") +
  labs(title = "PHÂN BỐ SỐ TRANG THEO LOẠI SP",
       x = NULL,
       y = "Số trang") +
  theme_minimal() +
  scale_fill_brewer(palette = "Pastel1") +
  theme(legend.position = "none")

# 2. Violin plot - Giữ nguyên
p2 <- ggplot(data, aes(x = LoaiSanPham, y = SoTrang, fill = LoaiSanPham)) +
  geom_violin(alpha = 0.5, trim = FALSE) +
  geom_jitter(width = 0.2, alpha = 0.3, size = 1.5) +
  labs(title = "MẬT ĐỘ PHÂN BỐ KÈM DATA POINTS",
       x = NULL,
       y = NULL) +
  theme_minimal() +
  scale_fill_brewer(palette = "Set2") +
  theme(legend.position = "none")

# 3. Barplot - Giữ nguyên
p3 <- ggplot(data, aes(x = LoaiSanPham, y = SoTrang, fill = LoaiSanPham)) +
  stat_summary(fun = mean, geom = "bar", position = "dodge") +
  stat_summary(fun.data = mean_sdl, 
               geom = "errorbar", 
               width = 0.2,
               color = "gray30") +
  labs(title = "GIÁ TRỊ TRUNG BÌNH ± ĐỘ LỆCH CHUẨN",
       x = NULL,
       y = "Số trang (Mean ± SD)") +
  theme_minimal() +
  scale_fill_brewer(palette = "Accent") +
  theme(legend.position = "none")

# 4. Sửa lại histogram phân nhóm
p4 <- ggplot(data, aes(x = SoTrang, fill = LoaiSanPham)) +
  geom_bar(stat = "count", alpha = 0.7, position = "dodge", color = "white") +
  facet_wrap(~LoaiSanPham, ncol = 1, scales = "free_y") +
  labs(title = "PHÂN BỐ TẦN SỐ CỦA SỐ TRANG",
       x = "Số trang",
       y = "Tần số") +
  theme_minimal() +
  scale_fill_brewer(palette = "Dark2") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) # Hiển thị nhiều giá trị hơn trên trục x

# Kết hợp 4 biểu đồ
combined_plot <- (p1 + p2) / (p3 + p4) + 
  plot_layout(heights = c(1, 1.5)) +
  plot_annotation(title = "THỐNG KÊ MÔ TẢ: LOẠI SP vs SỐ TRANG",
                  theme = theme(plot.title = element_text(size = 16, hjust = 0.5, face = "bold")))

# Xuất file
ggsave("LoaiSanPham_SoTrang.png", combined_plot, width = 12, height = 12, dpi = 300)

knitr::include_graphics("LoaiSanPham_SoTrang.png", error = FALSE)

data %>%
  filter(LoaiSanPham %in% c("Decal", "Biểu mẫu", "Catalog")) %>%
  count(LoaiSanPham, SoTrang) %>%
  group_by(LoaiSanPham) %>%
  gt() %>%
  tab_header(
    title = "BẢNG THỐNG KÊ TẦN SỐ CỦA SỐ TRANG",
    subtitle = "Phân theo Loại sản phẩm"
  ) %>%
  fmt_number(columns = n, decimals = 0) %>%
  cols_label(
    LoaiSanPham = "Loại sản phẩm",
    SoTrang = "Số trang",
    n = "Tần số"
  )

BẢNG THỐNG KÊ TẦN SỐ CỦA SỐ TRANG
Phân theo Loại sản phẩm
Số trang	Tần số
Biểu mẫu
80	38
81	78
82	84
83	79
84	76
85	83
86	88
87	67
88	87
89	87
90	98
91	79
92	66
93	71
94	71
95	90
96	88
97	87
98	84
99	82
100	48
Catalog
20	26
21	51
22	45
23	56
24	52
25	59
26	57
27	41
28	52
29	66
30	72
31	44
32	44
33	53
34	56
35	62
36	57
37	54
38	58
39	55
40	27
Decal
0	2,765

# Lọc Biểu mẫu và thống kê tần số
plot_freq <- data %>%
  filter(LoaiSanPham == "Biểu mẫu") %>%
  ggplot(aes(x = SoTrang)) +
  geom_bar(fill = "steelblue") +
  geom_text(stat = 'count', aes(label = after_stat(count)), 
            vjust = -0.5, size = 3) +
  labs(title = "TẦN SỐ SỐ TRANG BIỂU MẪU",
       x = "Số trang",
       y = "Tần số") +
  theme_minimal()

print(plot_freq)

# Lọc Catalog và thống kê tần số
plot_freq <- data %>%
  filter(LoaiSanPham == "Catalog") %>%
  ggplot(aes(x = SoTrang)) +
  geom_bar(fill = "steelblue") +
  geom_text(stat = 'count', aes(label = after_stat(count)), 
            vjust = -0.5, size = 3) +
  labs(title = "TẦN SỐ SỐ TRANG CATALOG",
       x = "Số trang",
       y = "Tần số") +
  theme_minimal()

print(plot_freq)

4.6 Thống kê mô tả các biến định lượng còn lại

summary(data$TyLeLoi)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.005001 0.009193 0.020083 0.017779 0.024927 0.029997

summary(data$SoLuongLoi)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.00    4.00   45.00   44.39   73.00  120.00

summary(data$DoanhThu)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  954800 1493100 1692000 1698486 1880200 3042000

summary(data$ChiPhiBienDong)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5920   12000   21340   21207   28712   46560

summary(profit_data$TongDoanhThu)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 17239300 23082300 25549300 25655086 28218600 34574100

summary(profit_data$TongChiPhiBienDong)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  160220  278811  318440  320324  360877  493791

summary(profit_data$LoiNhuan)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  2034995  7821102 10267809 10334762 12886244 19080309

5 Phân phối và diễn biến của Lợi nhuận theo thời gian

# 1. Phân phối lợi nhuận
ggplot(profit_data, aes(x = LoiNhuan)) +
  geom_histogram(bins = 30, fill = "steelblue", alpha = 0.7) +
  labs(title = "Phân phối Lợi nhuận hàng ngày",
       x = "Lợi nhuận (VND)", 
       y = "Tần số") +
  scale_x_continuous(labels = scales::comma) +  # Hiển thị số đầy đủ
  scale_y_continuous(labels = scales::comma) +  # Áp dụng cho cả trục y nếu cần
  theme_minimal()

# 2. Xu hướng thời gian
ggplot(profit_data, aes(x = Ngay, y = LoiNhuan)) +
  geom_line(color = "tomato") +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Diễn biến Lợi nhuận theo thời gian",
       x = "Ngày", 
       y = "Lợi nhuận (VND)") +
  scale_y_continuous(labels = scales::comma) +  # Hiển thị số đầy đủ
  theme_minimal()

6 Mô phỏng các vòng lặp của lợi nhuận

# Tạo hàm lấy mẫu ngẫu nhiên và tính lợi nhuận trung bình
simulate_profit <- function(n, profit_data) {
  means <- numeric(n)
  for (i in 1:n) {
    sampled_days <- sample(profit_data$LoiNhuan, size = n, replace = TRUE)
    means[i] <- mean(sampled_days)
  }
  return(means)
}

# Danh sách số vòng lặp
iterations <- c(10, 100, 1000, 10000)

# Thực hiện mô phỏng
simulations <- lapply(iterations, function(n) simulate_profit(n, profit_data))

# Chuẩn bị dữ liệu để vẽ biểu đồ
sim_df <- data.frame(
  Iteration = rep(iterations, times = sapply(simulations, length)),
  EstimatedProfit = unlist(simulations)
)

# Gán tên dễ hiểu
sim_df$Label <- paste0("n = ", sim_df$Iteration)

# Vẽ biểu đồ giống hình ảnh π
ggplot(sim_df, aes(x = Label, y = EstimatedProfit)) +
  geom_jitter(width = 0.1, height = 0, alpha = 0.4, color = "steelblue") +
  stat_summary(fun = mean, geom = "point", size = 5, color = "darkred") +
  labs(title = "Ước lượng lợi nhuận trung bình",
       x = "Số lần mô phỏng (n)",
       y = "Lợi nhuận trung bình/ngày (VND)") +
  scale_y_continuous(labels = scales::comma) +
  theme_minimal(base_size = 14)

# Hàm mô phỏng lợi nhuận trung bình
simulate_and_plot <- function(n, profit_data) {
  means <- numeric(n)
  for (i in 1:n) {
    sampled_days <- sample(profit_data$LoiNhuan, size = n, replace = TRUE)
    means[i] <- mean(sampled_days)
  }
  mean_est <- mean(means)
  
  # Tạo biểu đồ
  p <- ggplot(data.frame(x = means), aes(x = x)) +
    geom_histogram(bins = 30, fill = "#377eb8", alpha = 0.7, color = "black") +
    geom_vline(xintercept = mean_est, color = "red", linetype = "dashed") +
    labs(
      title = paste0(n, " Lần Mô Phỏng"),
      subtitle = paste0("Lợi nhuận TB ≈ ", round(mean_est, 0), " VND"),
      x = NULL,
      y = NULL
    ) +
    scale_x_continuous(labels = scales::comma) +
    theme_minimal(base_size = 13)
  return(p)
}

# Danh sách số vòng lặp
iterations <- c(10, 100, 1000, 10000)

# Tạo danh sách biểu đồ
plots <- lapply(iterations, function(n) simulate_and_plot(n, profit_data))

# Dùng patchwork để sắp xếp thành 2 hàng 3 cột
final_plot <- (plots[[1]] + plots[[2]]) / (plots[[3]] + plots[[4]])

# Hiển thị
ggsave("Mô phỏng Lợi nhuận.png", final_plot, width = 12, height = 6)

knitr::include_graphics("Mô phỏng Lợi nhuận.png", error = FALSE)

7 Các kịch bản của mô hình

# 1. HÀM MÔ PHỎNG LỢI NHUẬN
simulate_profit_mean <- function(profit_vector, n) {
  replicate(n, mean(sample(profit_vector, n, replace = TRUE)))
}

# 2. TẠO KỊCH BẢN MỚI (thay đổi dữ liệu)

create_modified_profit_data <- function(data, factor_error = 1, factor_cost = 1, factor_order = 1) {
  df <- data
  
  # Tăng/giảm tỷ lệ lỗi
  df$TyLeLoi <- pmin(df$TyLeLoi * factor_error, 1)  # Không vượt 100%
  df$SoLuongLoi <- round(df$SoLuongDatHang * df$TyLeLoi)
  
  # Tăng/giảm giá vốn
  df$GiaVon <- df$GiaVon * factor_cost
  
  # Tính lại chi phí biến động
  df$ChiPhiBienDong <- mapply(function(type, error_q, pages, cost) {
    if (type == "Decal") error_q * cost else error_q * pages * cost
  }, df$LoaiSanPham, df$SoLuongLoi, df$SoTrang, df$GiaVon)
  
  # Tăng/giảm số đơn hàng/ngày
  daily_orders <- df %>%
    group_by(Ngay) %>%
    summarise(
      TongDoanhThu = sum(DoanhThu, na.rm = TRUE),
      TongChiPhiBienDong = sum(ChiPhiBienDong, na.rm = TRUE),
      ChiPhiCoDinh = 15000000,
      .groups = "drop"
    ) %>%
    mutate(
      LoiNhuan = (TongDoanhThu * factor_order) - ChiPhiCoDinh - (TongChiPhiBienDong * factor_order)
    )
  
  return(daily_orders$LoiNhuan)
}

# 3. TẠO KỊCH BẢN VÀ MÔ PHỎNG
scenarios <- list(
  "Bình thường" = c(1, 1, 1),
  "Tăng lỗi" = c(1.5, 1, 1),
  "Giảm lỗi" = c(0.5, 1, 1),
  "Tăng giá vốn" = c(1, 1.2, 1),
  "Giảm giá vốn" = c(1, 0.8, 1),
  "Tăng đơn hàng" = c(1, 1, 1.3),
  "Giảm đơn hàng" = c(1, 1, 0.7)
)

iterations <- c(10, 100, 1000, 10000)

# Lưu kết quả mô phỏng
result_list <- list()

for (scenario_name in names(scenarios)) {
  factors <- scenarios[[scenario_name]]
  profit_vec <- create_modified_profit_data(data, 
                                            factor_error = factors[1], 
                                            factor_cost = factors[2], 
                                            factor_order = factors[3])
  
  for (n in iterations) {
    profits <- simulate_profit_mean(profit_vec, n)
    result_list[[paste0(scenario_name, "_n", n)]] <- data.frame(
      Scenario = scenario_name,
      Iteration = n,
      EstimatedProfit = profits
    )
  }
}

# Gộp dữ liệu
all_simulations <- bind_rows(result_list)

# 4. VẼ BIỂU ĐỒ
aa <- ggplot(all_simulations, aes(x = factor(Iteration), y = EstimatedProfit, fill = Scenario)) +
  geom_boxplot(alpha = 0.6, outlier.shape = NA) +
  facet_wrap(~Scenario, scales = "free_y") +
  labs(
    title = "Ảnh hưởng của các yếu tố đến lợi nhuận trung bình/ngày",
    x = "Số vòng lặp mô phỏng (n)",
    y = "Lợi nhuận trung bình/ngày (VND)"
  ) +
  theme_minimal(base_size = 13) +       
  theme(legend.position = "none")

print(aa)

ggsave("Ảnh hưởng LN.png", aa, width = 12, height = 6)

# Tính lợi nhuận trung bình/ngày theo từng kịch bản và từng số vòng lặp
summary_table <- all_simulations %>%
  group_by(Scenario, Iteration) %>%
  summarise(MeanProfit = mean(EstimatedProfit), .groups = "drop") %>%
  arrange(Iteration)

# Hiển thị
datatable(summary_table)

8 Đánh giá hiệu quả mô phỏng

evaluate_simulation <- function(sim_values, confidence = 0.95) {
  n <- length(sim_values)
  mean_val <- mean(sim_values)
  sd_val <- sd(sim_values)
  sem <- sd_val / sqrt(n)
  z <- qnorm((1 + confidence)/2)  # z = 1.96 với CI 95%
  ci_lower <- mean_val - z * sem
  ci_upper <- mean_val + z * sem
  cv <- sd_val / mean_val
  
  var_5 <- quantile(sim_values, probs = 0.05)
  es_5 <- mean(sim_values[sim_values <= var_5])
  
  return(data.frame(
    N = n,
    Mean = mean_val,
    SD = sd_val,
    SEM = sem,
    CI_Lower = ci_lower,
    CI_Upper = ci_upper,
    CV = cv,
    VaR_5 = var_5,
    ES_5 = es_5
  ))
}

# Áp dụng đánh giá cho từng kịch bản + n
evaluation_results <- all_simulations %>%
  group_by(Scenario, Iteration) %>%
  summarise(
    Eval = list(evaluate_simulation(EstimatedProfit)),
    .groups = "drop"
  ) %>%
  unnest(Eval) 

# Hiển thị bảng kết quả
datatable(evaluation_results)