用pymc实施MCMC

德国坦克问题的简单解答

德国坦克问题是一个经典的统计估计问题，不清楚的同学可以自行google。最简单的思路是将坦克总数X看做是一个X面的骰子。我们对骰子的面数不确定，但可以根据出现的点数来对X的取值分布做推断。推断方法就是贝叶斯方法。

假设我们对坦克总数有[500, 1000, 1500, 2000, 3000]五种可能的估计，先验估计为均匀分布，之后我们观察到战场上坦克编号出现有[4, 180, 75, 1007, 1003, 1500]

尝试估计这五种可能的后验概率分布如下：
500： 0.0
1000： 0.0
1500： 0.837
2000： 0.149
3000： 0.013

python 代码如下：

python贝叶斯文分类识别垃圾短信

Python贝叶斯文本分类识别垃圾短信

1、读取数据，type表示短信类别，text是短信内容

In [17]:

%pylab inline
import pandas as pd
import numpy as np
df = pd.read_csv('sms_spam.csv')
df.head()

Populating the interactive namespace from numpy and matplotlib

Out[17]:

	type	text
0	ham	Hope you are having a good week. Just checking in
1	ham	K..give back my thanks.
2	ham	Am also doing in cbe only. But have to pay.
3	spam	complimentary 4 STAR Ibiza Holiday or £10,000 ...
4	spam	okmail: Dear Dave this is your final notice to...

2、使用sklearn包转换文本为结构化数据，将矩阵分切为训练集和检验集

CountVectorizer负责将文档转为文档词频矩阵，重要的参数有如下几个：

• ngram_range：ngrame频率范围，如果需要识别词组的话需要设置
• stop_words：停词列表
• token_pattern：分词的字符模式，默认空格
• max_df：词频上限，超过该值的词项不作为特征，即过滤常用词
• min_df：词频下限，低于该值的词项不作为特征
• max_features：只选择词频较高的几个作为特征

In [18]:

from sklearn.feature_extraction.text import CountVectorizer
vectorizer = CountVectorizer(ngram_range=(1,1),stop_words='english',lowercase=True,min_df=1)
X = vectorizer.fit_transform(df.text) 
y = (df.type == 'spam').values.astype(int)

TfidfVectorizer则可以计算tfidf值，而非仅仅文档词频矩阵

In [19]:

from sklearn.feature_extraction.text import TfidfVectorizer
vectorizer = TfidfVectorizer(ngram_range=(1,1),stop_words='english',lowercase=True,min_df=1)
X = vectorizer.fit_transform(df.text)

3、将数据切分为train和test

In [20]:

from sklearn.cross_validation import train_test_split
xtrain, xtest, ytrain, ytest = train_test_split(X, y)

4、使用贝叶斯分类器进行训练

• 重要的参数alpha用于设置平滑系数

In [21]:

from sklearn.naive_bayes import MultinomialNB
clf = MultinomialNB(alpha =1).fit(xtrain, ytrain)

5、观察分类效果

In [22]:

training_accuracy = clf.score(xtrain, ytrain)
test_accuracy = clf.score(xtest, ytest)
print "训练集准确率:  {:.2f}".format(training_accuracy)
print "检验集准确率:  {:.2f}".format(test_accuracy)  

训练集准确率:  0.98
检验集准确率:  0.97

6、使用CV选择最优参数，参数为0.2

In [23]:

from sklearn import svm, grid_search
nb = MultinomialNB()
parameters = {'alpha':np.linspace(0,10,101)}
clf = grid_search.GridSearchCV(nb, parameters)
clf.fit(X, y)

Out[23]:

GridSearchCV(cv=None,
       estimator=MultinomialNB(alpha=1.0, class_prior=None, fit_prior=True),
       fit_params={}, iid=True, loss_func=None, n_jobs=1,
       param_grid={'alpha': array([  0. ,   0.1, ...,   9.9,  10. ])},
       pre_dispatch='2*n_jobs', refit=True, score_func=None, scoring=None,
       verbose=0)

In [100]:

print "最佳参数:  {:.3f}".format(clf.best_params_['alpha'] )
print "最佳准确率:  {:.3f}".format(clf.best_score_)  

最佳参数:  0.200
最佳准确率:  0.984

In [62]:

accuracy = [t[1] for t in clf.grid_scores_]
para = [t[0]['alpha'] for t in clf.grid_scores_]

In [94]:

import matplotlib.pylab as plt
accuracy = [t[1] for t in clf.grid_scores_]
para = [t[0]['alpha'] for t in clf.grid_scores_]
plt.plot(para,accuracy,lw=3)

Out[94]:

[]

理解MCMC

在贝叶斯统计中，经常需要计算后验概率，概率计算就涉及到积分问题。一种解决方法是用解析式得到后验概率直接计算，另一种是利用统计模拟来计算近似值。 考虑一个简单问题，我们对一个硬币反复投掷，对于出现正面的概率theta先主观设定为一个均匀分布，然后实际投掷14次，得到11次正面，要根据这个信息data来更新后验概率。下面用统计模拟来计算。 因为实际的数据表现是正面出现次数较多，所以后验概率会做出相应调整。这个统计模拟的例子非常简单，容易实施。但如果涉及的参数个数很多时，计算量就会非常大。此时就可以尝试另一种统计模拟方法，即MCMC(Markov chain Monte Carlo)。我们先考虑一个思想实验：

一个岛国下面有7个岛，各岛人数不一样。旅行者已经在其中一个岛上，他要去各岛游历，准备在人口多的岛上游玩的时间长一些。但他并不知道具体的人口数，可以询问本岛和相邻岛的市长。他的旅行计划如下：
首先每天会扔一枚硬币，如果正面向上，就计划去左岛，如果反面向上，就计划去右岛。
然后比较当前所在岛人口数a和计划去的目的岛人口数b，若b小于a，则实际去的概率为b/a，若b大于a，则实际去的概率为1。所以实际去的概率归纳为min(1,b/a)。
就这样，旅行者每天不断的在不同岛间移动，最终他在各岛上呆的时间之比，就会收敛为各岛人口之比，这样就通过一个随机模拟得到了一个概率分布。

朴素贝叶斯分类与贝叶斯网络

朴素贝叶斯分类（Naive Bayes Classifier）是一种简单而容易理解的分类方法，看起来很Naive，但用起来却很有效。其原理就是贝叶斯定理，从数据中得到新的信息，然后对先验概率进行更新，从而得到后验概率。好比说我们判断一个人的品质好坏，对于陌生人我们对他的判断是五五开，如果说他做了一件好事，那么这个新的信息使我们判断他是好人的概率增加了。朴素贝叶斯分类的优势在于不怕噪声和无关变量，其Naive之处在于它假设各特征属性是无关的。而贝叶斯网络（Bayesian Network）则放宽了变量无关的假设，将贝叶斯原理和图论相结合，建立起一种基于概率推理的数学模型,对于解决复杂的不确定性和关联性问题有很强的优势。下面我们用mlbench包中的一个数据集来看看如何用这两种方法进行学习训练。

PimaIndiansDiabetes2是美国一个疾病研究机构所拥有的一个数据集，其中包括了9个变量，共有768个样本。响应变量即是对糖尿病的判断，它是一个二元变量。其它各解释变量是个体的若干特征，如年龄和其它医学指标，均为数值变量。其中有一些缺失值的存在，虽然朴素贝叶斯分类对于缺失值并不敏感，我们还是先将其进行插补，再进行建模以观察准确率。R语言中的e1071包中就有可以实施朴素贝叶斯分类的函数，但在本例我们使用klaR包中的NaiveBayes函数，因为该函数较之前者增加了两个功能，一个是可以输入先验概率，另一个是在正态分布基础上增加了核平滑密度函数。为了避免过度拟合，在训练时还要将数据分割进行多重检验，所以我们还使用了caret包的一些函数进行配合。

# 加载扩展包和数据
library(caret)
data(PimaIndiansDiabetes2,package='mlbench')
# 对缺失值使用装袋方法进行插补
preproc <- preProcess(PimaIndiansDiabetes2[-9],method="bagImpute")
data <- predict(preproc,PimaIndiansDiabetes2[-9])
data$Class <- PimaIndiansDiabetes2[,9]
# 使用朴素贝叶斯建模，这里使用了三次10折交叉检验得到30个结果
fitControl <- trainControl(method = "repeatedcv", number = 10, repeats = 3,returnResamp = "all")
model1 <- train(Class~., data=data,method='nb',trControl = fitControl,tuneGrid = data.frame(.fL=1,.usekernel=F))
# 观察30次检验结果，发现准确率在0.75左右
resampleHist(model1)
# 返回训练数据的混淆矩阵
pre <- predict(model1)
confusionMatrix(pre,data$Class)