Daily Income of Sex Workers in Hanoi and Saigon

R for Fun

Nguyen Chi Dung

Thông tin về giá của dịch vụ (và đánh giá của của người sử dụng) lấy từ https://gaigoivl.com/.

R Codes để cào dữ liệu ở đây.

Nếu không muốn mất thời gian cho việc cào dữ liệu thì lấy data ở đây.

library(magrittr)
library(stringr)
library(tidyverse)
rm(list = ls())

# Đọc dữ liệu: 

pho <- read.csv("D:/du_lieu_pho.csv")

# 133 quan sát đầu là thuộc khu vực HN, còn lại ở SG: 
pho %<>% mutate(city = c(rep("HN", 133), rep("SG", nrow(.) - 133)))
pho %<>% mutate_if(is.factor, as.character)

# Viết xử lí thông tin về giá: 

get_price <- function(x) {
  x %>% 
    str_replace_all("[^0-9]", "") %>% 
    as.numeric() %>% 
    return()
}

# Có thể thấy mức giá phổ biến là 400 - 500 cho một lần sử dụng dịch vụ: 

pho %<>% mutate(price = get_price(price))

theme_set(theme_minimal())
pho %>% 
  group_by(price) %>% 
  count() %>% 
  na.omit() %>% 
  ggplot(aes(price, n)) + 
  geom_col() + 
  labs(x = NULL, 
       y = "Price", 
       title = "Price Distribution of Sex Service", 
       caption = "Data Source: https://gaigoivl.com/")

# Theo từng thành phố: 

pho %>% 
  group_by(price, city) %>% 
  count() %>% 
  na.omit() %>% 
  ggplot(aes(price, n)) + 
  geom_col() + 
  facet_wrap(~ city, scales = "free") + 
  labs(x = NULL, 
       y = "Price", 
       title = "Price Distribution of Sex Service by City", 
       caption = "Data Source: https://gaigoivl.com/")

# Giá của dịch vụ sex trung bình ở SG cao hơn HN: 

pho %>% 
  filter(!is.na(price)) %>% 
  group_by(city) %>% 
  summarise_each(funs(mean, median, min, max, n()), price) %>% 
  mutate_if(is.numeric, function(x) {round(x, 0)}) %>% 
  knitr::kable(col.names = c("City", "Mean", "Median", "Min", "Max", "N"))

# Mô phỏng Monte Carlo cho thu nhập ngày của các em phò với
# giả định then chốt rằng tần suất phục vụ mỗi ngày của 
# các em là từ 2 đến 5 và phân phố này là Uniform: 


income_simu <- function(n) {
  price <- pho$price %>% na.omit()
  k <- length(price)
  income_si <- c()
  for (i in 1:n) {
    tan_suat_phuc_vu <- sample(2:5, size = k, replace = TRUE)
    income <- mean(price*tan_suat_phuc_vu)
    income_si <- c(income_si, income)
  }
  return(income_si)
}

# Mô phỏng với 50 ngàn lần: 

daily_income <- income_simu(50000)
daily_income <- data.frame(income = daily_income)

# Có dấu hiệu để tin rằng thu nhập theo ngày của các 
# em phò là phân phối chuẩn với trung bình khoảng 1.6 triệu
# và gần như chắc chắn rằng thu nhập ngày của các em phò
# nằm đâu đó trong khoảng từ 1493k đến 1720k. 

daily_income %>% 
  ggplot(aes(income)) + 
  geom_density(alpha = 0.3, fill = "blue", color = "blue") + 
  geom_histogram(aes(y = ..density..), fill = "red", color = "red", alpha = 0.3) + 
  geom_vline(xintercept = mean(daily_income$income), color = "green", size = 1.1) + 
  labs(x = NULL, y = NULL, 
       subtitle = "Based on Assumption that Daily Service Frequency is from 2 to 5 and Distribution is Uniform", 
       title = "Daily Income of Sex Workers based on Naive Monte Carlo Simulation", 
       caption = "Data Source: https://gaigoivl.com")

daily_income$income %>% summary()

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1488    1589    1607    1607    1626    1719