hw4
黄玺竹
2024-11-02
Problem 1
最大似然函数： \[\begin{align} l(\theta)=&log\prod_{i=1}^N P(x_i,y_i;\theta)\\ =&log\prod_{i=1}^N P(x_i|y_i;\theta)P(y_i;\theta)\\ =&log\prod_{i=1}^N \bigg(\prod_{j=1}^pP(x_i^{(j)}|y_i;\theta)\bigg)P(y_i;\theta)\\ =&\sum_{i=1}^N\bigg[\sum_{j=1}^plogP(x_i^{(j)}|y_i;\theta)+logP(y_i;\theta)\bigg] \end{align}\] 其中， \[P(Y=y_i)=\prod_{k=1}^KP(Y=c_k)^{I(y_i=c_k)}\\ P(X^{(j)}=x_i^{(j)}|Y=y_i)=\prod_{l=1}^{S_j}P(X^{(j)}=a_{jl}|Y=c_k)^{I(x_i^{(j)}=a_{jl},~y_i=c_k)}\]
所以最大似然函数为 \(l(\theta)=\sum_{i=1}^N\bigg[\sum_{j=1}^p\sum_{l=1}^{S_j}I(x_i^{(j)}=a_{jl},y_i=c_k)logP(X^{(j)}=a_{jl}|Y=c_k)+\sum_{k=1}^KI(y_i=c_k)logP(Y=c_k)\bigg]\) 需要满足以下约束条件： \[\sum_{k=1}^KP(Y=c_k)=1\\ \sum_{l=1}^{S_j}P(X^{(j)}=a_{jl}|Y=c_k)=1\quad for ~every~k=1,\cdots,K. \]
可以得到\(P\left(Y=c_k\right)\)的拉格朗日函数为： \[ L(\lambda) = \sum_{i=1}^N \sum_{k=1}^K I\left(y_i=c_k\right) \log P\left(Y=c_k\right) + \lambda \left( \sum_{k=1}^{K}P( Y = c_k)- 1 \right) \]
\[ \frac{\partial L(\lambda)}{\partial P\left(Y=c_k\right)} = \sum_{i=1}^{N}\frac{I\left(y_i=c_k\right)}{P\left(Y=c_k\right)} +\lambda = 0 \] \[\frac{\partial L(\lambda)}{\partial \lambda} =\sum_{k=1}^KP(Y=c_k)-1=0\] 根据上面两式可得： \[\sum_{i=1}^NI(y_i=c_k)+\lambda P(Y=c_k)=0\\ \sum_{k=1}^K\sum_{i=1}^NI(y_i=c_k)+\lambda \sum_{k=1}^KP(Y=c_k)=0\] 因此 \[\lambda=-\sum_{k=1}^K\sum_{i=1}^NI(y_i=c_k).\] \[P(Y=c_k)=-\frac{\sum_{i=1}^N I(y_i=c_k)}{\lambda}=\frac{\sum_{i=1}^N I(y_i=c_k)}{\sum_{k=1}^K\sum_{i=1}^NI(y_i=c_k)}\] 显然\(\sum_{k=1}^K\sum_{i=1}^NI(y_i=c_k)=N\).
因此\[P(Y=c_k)=\frac{\sum_{i=1}^N I(y_i=c_k)}{N}\] 同理可得到 \(P(X^{(j)}=a_{j l} \mid Y=c_k)\) 的拉格朗日函数为： \[ L(\mu) = \sum_{i=1}^{N}\sum_{j=1}^n \sum_{l=1}^{S_j} I\left(x_i^{(j)}=a_{j l}, y_i=c_k\right) \log P\left(X^{(j)}=a_{j l} \mid Y=c_k\right) + \mu\left( \sum_{l=1}^{S_j} P\left(X^{(j)}=a_{j l} \mid Y=c_k\right) - 1 \right) \] \[ \frac{\partial L(\lambda)}{\partial P\left(X^{(j)}=a_{j l} \mid Y=c_k\right)} = \sum_{i=1}^{N} \frac{I\left(x_i^{(j)}=a_{j l}, y_i=c_k\right)}{P\left(X^{(j)}=a_{j l} \mid Y=c_k\right)} + \mu = 0\] \[\frac{\partial L(\mu)}{\partial \mu} =\sum_{l=1}^{S_j}P(X^{(j)}=a_{j l} \mid Y=c_k)-1=0\] 根据上面两式可得： \[\mu=-\sum_{l=1}^{S_j}\sum_{i=1}^NI(x_i^{(j)}=a_{jl},y_i=c_k)=-\sum_{i=1}^NI(y_i=c_k)\] \[P(X^{(j)}=a_{j l} \mid Y=c_k)=-\frac{\sum_{i=1}^NI(x_i^{(j)}=a_{jl},y_i=c_k)}{\mu}=\frac{\sum_{i=1}^NI(x_i^{(j)}=a_{jl},y_i=c_k)}{\sum_{i=1}^NI(y_i=c_k)}\]
Problem 2
MLE
先验概率：\(P(Y=c_k)=\frac{\sum_{i=1}^NI(y_i=c_k)}{N}\)
条件概率：\(P(X^{(j)}=a_{j l} \mid Y=c_k)=\frac{\sum_{i=1}^NI(x_i^{(j)}=a_{jl},y_i=c_k)}{\sum_{i=1}^NI(y_i=c_k)}\)
所以对于\(x=(2,~M)^T\)， \[P(Y=1)=\frac{2}{3},\quad P(Y=-1)=\frac{1}{3}\\ P(X^{(1)}=2\mid Y=1)=\frac{1}{5},\quad P(X^{(1)}=2\mid Y=-1)=\frac{3}{5}\\ P(X^{(2)}=M\mid Y=1)=\frac{2}{5},\quad P(X^{(2)}=M\mid Y=-1)=\frac{2}{5} \] \[\begin{align} y=&\arg\max_{c_k}P(Y=c_k)\prod_{j=1}^2P(X^{(j)}=x^{(j)}\mid Y=c_k)\\ =&\arg\max_{c_k}\quad \{\frac{4}{75},~\frac{6}{75}\}\\ =& -1 \end{align}\]
贝叶斯估计
先验概率：\(P(Y=c_k)=\frac{\sum_{i=1}^NI(y_i=c_k)+\lambda}{N+\lambda K}\)
条件概率：\(P(X^{(j)}=a_{j l} \mid Y=c_k)=\frac{\sum_{i=1}^NI(x_i^{(j)}=a_{jl},y_i=c_k)+\lambda}{\sum_{i=1}^NI(y_i=c_k)+S_j \lambda}\)
所以对于\(x=(2,~M)^T\)， \[P(Y=1)=\frac{10+1}{15+2}=\frac{11}{17},\quad P(Y=-1)=\frac{5+1}{15+2}=\frac{6}{17}\\ P(X^{(1)}=2\mid Y=1)=\frac{2+1}{10+3}=\frac{3}{13},\quad P(X^{(1)}=2\mid Y=-1)=\frac{3+1}{5+3}=\frac{1}{2}\\ P(X^{(2)}=M\mid Y=1)=\frac{4+1}{10+3}=\frac{5}{13},\quad P(X^{(2)}=M\mid Y=-1)=\frac{2+1}{5+3}=\frac{3}{8} \] \[\begin{align} y=&\arg\max_{c_k}P(Y=c_k)\prod_{j=1}^2P(X^{(j)}=x^{(j)}\mid Y=c_k)\\ =&\arg\max_{c_k}\quad \{\frac{11}{17}×\frac{3}{13}×\frac{5}{13},~\frac{6}{17}×\frac{1}{2}×\frac{3}{8}\}\\ =&\arg\max_{c_k}\quad\{\frac{165}{2873},~\frac{18}{272}\}\\ =& -1 \end{align}\]
Problem 3
(1)
读入数据，命名为catalog。使用head()函数查看数据的基本形式。
# 设置工作路径
setwd("D:/File/大三上/统计学习概论/homework/homework_4")
# 读入catalog数据
catalog = read.csv('catalogs.csv', encoding = "utf-8")
# 呈现前六行数据
head(catalog)## name first
## 1 双汇 猪舌 700g/袋 生鲜
## 2 （满38减10）黑猪腊肉 湘西腊肉 即食小吃食猪肉干肉脯湖南张家界土家 香辣味 生鲜
## 3 精气神 猪肋排段（排骨） 400g/袋 山黑猪 黑猪肉 林间散养 生鲜
## 4 双汇 猪五花肉片 300g/袋 整肉原切 火锅食材 生鲜
## 5 农家散养新鲜现杀黑土猪肉 纯瘦肉 500g 生鲜
## 6 精气神 猪腿肉 400g/袋 山黑猪 黑猪肉 林间散养（2件起售） 生鲜
## second third
## 1 猪牛羊肉 猪肉
## 2 猪牛羊肉 猪肉
## 3 猪牛羊肉 猪肉
## 4 猪牛羊肉 猪肉
## 5 猪牛羊肉 猪肉
## 6 猪牛羊肉 猪肉针对二级分类“粮油调味”下的各三级分类分布情况绘制饼状图
##加plyr包
library(plyr)
# 读取粮油调味类产品数据并按照产品数量排序
df1 = ddply(catalog[catalog$second=="粮油调味",], .(third), nrow)
df1 = df1[order(df1$V1, decreasing = TRUE),]
# 建立数据标签
label = as.vector(df1$third)
label = paste(label, "(", round(df1$V1/sum(df1$V1)*100, 2), "%)", sep = "")
# R自带函数饼图绘制
pie(df1$V1,labels = label,main = "粮油调味类产品下三级分类饼图")
解读如下：
1.【粮油调味】下共有6个三级品类，其中调味品、方便食品、南北干货的总占比约75%，是该粮油调味品类下占比较高的三级品类；
2.除此之外，米面杂粮和食用油占比相当（11%左右），而有机食品占比最少，只有2.26%。
（2）
library(jiebaR)## 载入需要的程序包：jiebaRDlibrary(stringr)
##初始化分词工具
seg = worker(bylines = T)
##读入自定义词典，并加入到分词引擎中，注意到自定义词典中的名称也需要分词
new_words <- read.table('userdict.dat', header = F, stringsAsFactors = F)
new_user_word(seg, new_words$V1)## [1] TRUE# 数据文本化并进行分词
catalog$name = as.character(catalog$name)
cata_seg = segment(catalog$name, seg)
grain_and_oil_seg = cata_seg[which(catalog$second == "粮油调味")]
# 计算不同商品名称包含的词数
word_lengths = sapply(grain_and_oil_seg, length)
hist(word_lengths)
summary(word_lengths)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 12.00 15.00 15.17 18.00 35.00library(wordcloud2)
# 建立分词后词频数据框
grain_and_oil_seg_df = data.frame(table(unlist(grain_and_oil_seg)))
# 对词频进行降序排列，找到高频词
grain_and_oil_seg_df = grain_and_oil_seg_df[order(grain_and_oil_seg_df$Freq,decreasing = T),]
head(grain_and_oil_seg_df)## Var1 Freq
## 590 g 5892
## 2177 袋 1488
## 367 5 1400
## 163 2 1152
## 2437 调料 1110
## 374 500 1070wordcloud2(grain_and_oil_seg_df[1:100,])1.粮油调料类商品的分词数量的平均值为15.17个，中位数为15，最少2个，最多35个，根据直方分布图来看，商品词数呈现右偏分布，绝大多数的商品名词数都在25词以内。
2.“g”、“袋”、“L”、“kg”、“ml”、“瓶”、“装”等词汇以及一些数字的出现频率较高，这些词汇通常与商品的包装、规格和单位有关，反映了粮油调味商品在描述时常用的单位和包装方式。但由于其通用性，不适合作为分类使用的特征。
3.“食用油”、“方便面”、“特产”、“调味料”、“干货”、“酱”、“下饭菜”等词汇的出现，表明了“粮油调味”类别中包含的商品类型，这些商品是粮油调味类别中的热门商品。
（3）
将分词后的数据集按照7:3比例随机拆分为训练集和测试集。对分词结果进行预处理，预处理包括去掉数字字母组合，然后去掉低频词（包含该词的商品数量<10）、高频词（在多于75%的商品名称中出现）.最后，保留高频词中的top1000的词语，作为分类器的特征，并建立文档-词频矩阵。
library(text2vec)
# 数据集拆分
set.seed(1000)
id <- sample(1:nrow(catalog), 0.7 * nrow(catalog), F)
train <- catalog[id, ]
test <- catalog[-id, ]
train$id <- as.character(1:nrow(train))
test$id <- as.character(1:nrow(test))
# 分词结果相应进行拆分
seg_train <- cata_seg[id]
seg_test <- cata_seg[-id]
# 去除停用词，筛选特征词，生成语料文件
# 找到分词结果中的纯数字或字母，存入stop_words
stop_words <- unique(unlist(seg_train))
stop_words <- stop_words[grepl('^([a-z]|[A-Z])+$', stop_words)|grepl('^[0-9]+(\\.)?[0-9]*(-)?[0-9]*(\\.)?[0-9]*$', stop_words)]
# stop_words中可能存在一些包含在自定义词典中的专有名词，去掉
stop_words <- setdiff(stop_words, new_words)
# 设置分词迭代器
it_train <- itoken(seg_train, ids = train$id)
it_test <- itoken(seg_test, ids = test$id)
vocab <- create_vocabulary(it_train)
# 修剪低频词、高频词
pruned_vocab <- prune_vocabulary(vocab, doc_count_min = 10, doc_proportion_max = 0.75)
pruned_vocab <- pruned_vocab[!(pruned_vocab$term %in% stop_words), ]
# 保留词频前1000的词语
pruned_vocab <- pruned_vocab[order(pruned_vocab$term_count, decreasing = T), ]
pruned_vocab <- pruned_vocab[1:1000, ]
# 初步形成语料文件
vectorizer <- vocab_vectorizer(pruned_vocab)
# 根据语料文件生成文档-词频矩阵
# 文档-词频矩阵
dtm_train <- create_dtm(it_train, vectorizer)
dtm_test <- create_dtm(it_test, vectorizer)
dtm_train <- as.matrix(dtm_train)
dtm_test <- as.matrix(dtm_test)
# 转化为0-1矩阵
dtm_train[dtm_train > 1] <- 1
dtm_test[dtm_test > 1] <- 1
dtm_train <- apply(dtm_train, 2, as.factor)
dtm_test <- apply(dtm_test, 2, as.factor)
dim(dtm_train)## [1] 48645 1000dim(dtm_test)## [1] 20849 1000以文档词频矩阵为自变量，以每个商品名称对应的三级分类为因变量，对训练集建立朴素贝叶斯分类器。
library(e1071)
train_grain <- train[train$second == '粮油调味', ]
test_grain <- test[test$second == '粮油调味', ]
dtm_train_grain <- dtm_train[train$second == '粮油调味', ]
dtm_test_grain <- dtm_test[test$second == '粮油调味', ]
# 根据训练集数据，建立朴素贝叶斯分类器
nb <- naiveBayes(x = dtm_train_grain, y = as.factor(train_grain$third), laplace = 0)
# 在测试集上做预测
pred <- predict(nb, dtm_test_grain, type = 'class')
# 计算准确率
accuracy <- sum(as.character(test_grain$third) == as.character(pred)) / nrow(test_grain)
print(paste0('三级品类：', accuracy))## [1] "三级品类：0.905863708399366"