第八课、决策树

# 决策树(Decision Tree）是在已知各种情况发生概率的基础上，
# 通过构成决策树来求取净现值的期望值大于等于零的概率，评价项目风险，
# 判断其可行性的决策分析方法，是直观运用概率分析的一种图解法

# ### Demo: 选择给定功能的理想切点

# 车辆数据
import pandas as pd
url = 'https://raw.githubusercontent.com/justmarkham/DAT8/master/data/vehicles_train.csv';;
train = pd.read_csv(url)

# 在分割任何东西之前，只需预测整个数据集的平均值即可
train['prediction'] = train.price.mean()
train

#计算这些预测的RMSE（它的计算方法是先平方、再平均、然后开方）说明样本的离散程度。
from sklearn import metrics
import numpy as np
np.sqrt(metrics.mean_squared_error(train.price, train.prediction))

# 定义一个函数来计算一个给定的英里数的RMSE
def mileage_split(miles):
    lower_mileage_price = train[train.miles < miles].price.mean()
    higher_mileage_price = train[train.miles >= miles].price.mean()
    train['prediction'] = np.where(train.miles < miles, lower_mileage_price, higher_mileage_price)
    return np.sqrt(metrics.mean_squared_error(train.price, train.prediction))

# 计算分裂在英里数<50000的树的RMSE
print 'RMSE:', mileage_split(50000)
train

# 计算分裂在英里数<100000的树的RMSE
print 'RMSE:', mileage_split(100000)
train

# 检查所有可能的里程分割
mileage_range = range(train.miles.min(), train.miles.max(), 1000)
RMSE = [mileage_split(miles) for miles in mileage_range]

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (6, 4)
plt.rcParams['font.size'] = 14

# 绘制里程分界点（x轴）与RMSE（y轴）
plt.plot(mileage_range, RMSE)
plt.xlabel('Mileage cutpoint')
plt.ylabel('RMSE (lower is better)')

#**总结：**在每次分割之前，对每个特征重复该过程，并选择产生最低MSE的特征和分割点。
###在scikit-learn中构建一个回归树

#编码为0，卡车为1
train['vtype'] = train.vtype.map({'car':0, 'truck':1})

# define X and y
feature_cols = ['year', 'miles', 'doors', 'vtype']
X = train[feature_cols]
y = train.price

# 实例化一个DecisionTreeRegressor（random_state = 1）
from sklearn.tree import DecisionTreeRegressor
treereg = DecisionTreeRegressor(random_state=1)
treereg

# 使用leave-one-out交叉验证（LOOCV）来估计此模型的RMSE
from sklearn.cross_validation import cross_val_score
scores = cross_val_score(treereg, X, y, cv=14, scoring='neg_mean_squared_error')
np.mean(np.sqrt(-scores))

###调整回归树

#我们试着通过调整** max_depth **参数来减少RMSE：

#尝试不同的值一个接一个
treereg = DecisionTreeRegressor(max_depth=1, random_state=1)
scores = cross_val_score(treereg, X, y, cv=14, scoring='neg_mean_squared_error')
np.mean(np.sqrt(-scores))

#或者，我们可以写一个循环来尝试一系列值：
#要尝试的值列表
max_depth_range = range(1, 8)

# 列表来存储每个max_depth值的平均RMSE
RMSE_scores = []

# 使用每个值为max_depth的LOOCV
for depth in max_depth_range:
    treereg = DecisionTreeRegressor(max_depth=depth, random_state=1)
    MSE_scores = cross_val_score(treereg, X, y, cv=14, scoring='neg_mean_squared_error')
    RMSE_scores.append(np.mean(np.sqrt(-MSE_scores)))

# plot max_depth（x轴）对RMSE（y轴）
plt.plot(max_depth_range, RMSE_scores)
plt.xlabel('max_depth')
plt.ylabel('RMSE (lower is better)')

# max_depth = 3是最好的，所以适合使用该参数的树
treereg = DecisionTreeRegressor(max_depth=3, random_state=1)
treereg.fit(X, y)

# 每个特征的“基尼重要性”：该特征带来的（归一化）总误差减少
pd.DataFrame({'feature':feature_cols, 'importance':treereg.feature_importances_})

###创建一个树形图

#创建一个Graphviz文件
from sklearn.tree import export_graphviz
export_graphviz(treereg, out_file='tree_vehicles.dot', feature_names=feature_cols)

# ## 预测测试数据

# 阅读测试数据
url = 'https://raw.githubusercontent.com/justmarkham/DAT8/master/data/vehicles_test.csv'
test = pd.read_csv(url)
test['vtype'] = test.vtype.map({'car':0, 'truck':1})
test

#**问题：**使用上面的树形图，模型对每个观测值做出什么样的预测？
#使用拟合的模型来预测测试数据
X_test = test[feature_cols]
y_test = test.price
y_pred = treereg.predict(X_test)
y_pred

# 计算RMSE
np.sqrt(metrics.mean_squared_error(y_test, y_pred))

# 为你自己的树计算RMSE！
y_test = [3000, 6000, 12000]
y_pred = [0, 0, 0]
from sklearn import metrics
np.sqrt(metrics.mean_squared_error(y_test, y_pred))

###第2部分：分类树

# 读取数据
url = 'https://raw.githubusercontent.com/justmarkham/DAT8/master/data/titanic.csv';;
titanic = pd.read_csv(url)

# 将女性编码为0，男性编码为1
titanic['Sex'] = titanic.Sex.map({'female':0, 'male':1})

# 用年龄中位数填写年龄缺失值
titanic.Age.fillna(titanic.Age.median(), inplace=True)

# 创建一个虚拟变量的数据框
embarked_dummies = pd.get_dummies(titanic.Embarked, prefix='Embarked')
embarked_dummies.drop(embarked_dummies.columns[0], axis=1, inplace=True)

# 连接原始DataFrame和虚拟DataFrame
titanic = pd.concat([titanic, embarked_dummies], axis=1)

# 打印更新的DataFrame
titanic.head()

# - **存活：** 0 =死亡，1 =存活（反应变量）
# - ** Pclass：** 1 =一等，2 =二等，3 =三等
# - 如果树在这个特性上分裂会发生什么？
# - **性别：** 0 =女性，1 =男性
# - **年龄：**数值
# - **进入：** C或Q或S

# define X and y
feature_cols = ['Pclass', 'Sex', 'Age', 'Embarked_Q', 'Embarked_S']
X = titanic[feature_cols]
y = titanic.Survived

# 在所有数据上都适合一个max_depth = 3的分类树
from sklearn.tree import DecisionTreeClassifier
treeclf = DecisionTreeClassifier(max_depth=3, random_state=1)
treeclf.fit(X, y)

#创建一个Graphviz文件
export_graphviz(treeclf, out_file='tree_titanic.dot', feature_names=feature_cols)

# 计算特征的重要性
pd.DataFrame({'feature':feature_cols, 'importance':treeclf.feature_importances_})

##第3部分：将决策树与其他模型进行比较
# **决策树的优势：**
#  - 可以用于回归或分类
#  - 可以用图形显示
#  - 高度可解释的
#  - 可以指定为一系列规则，并且比其他模型更接近于人类决策
#  - 预测很快
#  - 功能不需要缩放
#  - 自动学习功能交互
#  - 倾向于忽略不相关的功能
#  - 非参数（如果特征和响应之间的关系非常非线性，将优于线性模型）

# ＃**决策树的缺点：**
# ＃
#  - 表现（通常）与最好的监督式学习方法没有竞争力
#  - 可以轻松过度训练数据（需要调整）
#  - 数据的小变化可能导致完全不同的树（高变化）
# 递归二进制分裂使得“局部最优”决策可能不会导致全局最优树
#  - 如果课程非常不平衡，则不会趋于良好
#  - 对于非常小的数据集不适合

第九课、集合

# 为什么我们学习集合？
#  - 提高机器学习模型预测性能的非常流行的方法
#  - 为理解更复杂的模型奠定了基础

# 第一部分：简介
import numpy as np

# 为可再现性设置种子
np.random.seed(1234)

#为每个模型生成1000个随机数（0到1之间），代表1000个观察值
mod1 = np.random.rand(1000)
mod2 = np.random.rand(1000)
mod3 = np.random.rand(1000)
mod4 = np.random.rand(1000)
mod5 = np.random.rand(1000)

# 如果随机数至少为0.3，则每个模型独立预测1（“正确响应”）
preds1 = np.where(mod1 > 0.3, 1, 0)
preds2 = np.where(mod2 > 0.3, 1, 0)
preds3 = np.where(mod3 > 0.3, 1, 0)
preds4 = np.where(mod4 > 0.3, 1, 0)
preds5 = np.where(mod5 > 0.3, 1, 0)

# 打印每个模型的前20个预测
print preds1[:20]
print preds2[:20]
print preds3[:20]
print preds4[:20]
print preds5[:20]

# 平均预测，然后轮到0或1
ensemble_preds = np.round((preds1 + preds2 + preds3 + preds4 + preds5)/5.0).astype(int)

#打印合奏的前20个预测
print ensemble_preds[:20]

#每个模型有多精确？
print preds1.mean()
print preds2.mean()
print preds3.mean()
print preds4.mean()
print preds5.mean()

#合奏有多精确？
print ensemble_preds.mean()

###什么是合成？
#**集合学习（或“集合”）**是将多个预测模型组合在一起以生成比任何单个模型更精确的组合模型的过程。

# 第2部分：手动合成
# 什么使得一个好的手工合奏？
#  - 不同类型的**模型**
#  - **功能的不同组合
#  - 不同的**调整参数**

# ＃＃第三部分：装袋
#
# **决策树**的主要弱点是它们往往不具有最好的预测准确性。 这部分是由于**高方差**，
# 这意味着训练数据中的不同分裂可能导致非常不同的树。
# ＃** Bagging **是减少机器学习方法方差的通用程序，但对决策树特别有用。
#  Bagging是** bootstrap聚合**的缩写，意思是bootstrap样本的聚合。

# 为可再现性设置种子
np.random.seed(1)

# 创建一个1到20的数组
nums = np.arange(1, 21)
print nums

# 样本数组替换20次
print np.random.choice(a=nums, size=20, replace=True)

# ＃**装袋工作（对于决策树）如何工作？**
# ＃1.使用训练数据中的B bootstrap样本生成B树。
# ＃2.训练每个树的自举样本并做出预测。
# ＃3.结合预测：
# ＃ - 平均**回归树的预测**
# ＃ - 为**分类树**投票

# ＃注意：
#  - **每个引导样本**应该与原始训练集大小相同。
#  - ** B **应该是一个足够大的值，这个错误似乎已经“稳定”了。
#  - 树木长得很深，所以它们的偏差很小/偏差很大。
# Bagging通过**降低方差来提高预测精度**，类似于交叉验证如何通过分割多次平均结果来减少与列车/测试分割相关联的方差（用于估计样本外误差）。

# ##手工实现袋装决策树（B = 10）

# 阅读并准备车辆训练数据
import pandas as pd
url = 'https://raw.githubusercontent.com/justmarkham/DAT8/master/data/vehicles_train.csv';;
train = pd.read_csv(url)
train['vtype'] = train.vtype.map({'car':0, 'truck':1})
train

#为可再现性设置种子
np.random.seed(123)

# 创建10个引导样本（将用于从DataFrame中选择行）
samples = [np.random.choice(a=14, size=14, replace=True) for _ in range(1, 11)]
samples

# 显示第一个决策树的行
train.iloc[samples[0], :]

#阅读并准备车辆测试数据
url = 'https://raw.githubusercontent.com/justmarkham/DAT8/master/data/vehicles_test.csv';;
test = pd.read_csv(url)
test['vtype'] = test.vtype.map({'car':0, 'truck':1})
test

from sklearn.tree import DecisionTreeRegressor

#长出每棵树
treereg = DecisionTreeRegressor(max_depth=None, random_state=123)

# 列表中存储每棵树的预测价格
predictions = []

# 定义测试数据
X_test = test.iloc[:, 1:]
y_test = test.iloc[:, 0]

# 为每个bootstrap样本生成一棵树，并对测试数据进行预测
for sample in samples:
    X_train = train.iloc[sample, 1:]
    y_train = train.iloc[sample, 0]
    treereg.fit(X_train, y_train)
    y_pred = treereg.predict(X_test)
    predictions.append(y_pred)

# 将预测从列表转换为NumPy数组
predictions = np.array(predictions)
predictions

# 平均预测
np.mean(predictions, axis=0)

# 计算RMSE
from sklearn import metrics
y_pred = np.mean(predictions, axis=0)
np.sqrt(metrics.mean_squared_error(y_test, y_pred))

# ## 用scikit-learn（B = 500）

# 定义培训和测试集
X_train = train.iloc[:, 1:]
y_train = train.iloc[:, 0]
X_test = test.iloc[:, 1:]
y_test = test.iloc[:, 0]

# 指示BaggingRegressor使用DecisionTreeRegressor作为“基本估计器”
from sklearn.ensemble import BaggingRegressor
bagreg = BaggingRegressor(DecisionTreeRegressor(), n_estimators=500, bootstrap=True, oob_score=True, random_state=1)

# 适合和预测
bagreg.fit(X_train, y_train)
y_pred = bagreg.predict(X_test)
y_pred

# 计算RMSE
np.sqrt(metrics.mean_squared_error(y_test, y_pred))

# ##估计样本外错误
# 对于袋装模型，可以不使用**火车/测试分割**或**交叉验证**来估计样本外错误！
# 平均而言，每个袋装树使用大约三分之二**的观察值。 对于每一棵树，其余的**观测**被称为“外包”观测。

# 显示第一个bootstrap示例
samples[0]

# 显示每个样品的“袋内”观察结果
for sample in samples:
    print set(sample)

# 显示每个样本的“外包”观察结果
for sample in samples:
    print sorted(set(range(14)) - set(sample))

# 如何计算**“外包错误”：**
# 1.对于训练数据中的每一个观测值，只使用**来预测其观测值在树袋外的树的响应值。 平均那些预测（回归）或投票（分类）。
# 2.将所有的预测与实际响应值进行比较，以便计算出袋外误差。
# 当B足够大时，**出包错误**是对样本外错误**的准确估计。

# 计算B = 500的袋外R平方得分（不幸的是，不是MSE）
bagreg.oob_score_

# ##估计特征重要性
# 套袋增加**预测的准确性**，但减少**模型的解释性**，因为不再可能将树形象化以理解每个特征的重要性。

# ＃第5部分：构建和调整决策树和随机森林
#  - 每个观察代表一个球员
#  - **进球：**预测球员薪水

###准备数据
url = 'https://raw.githubusercontent.com/justmarkham/DAT8/master/data/hitters.csv';;
hitters = pd.read_csv(url)

# 删除缺少值的行
hitters.dropna(inplace=True)

hitters.head()

# 将分类变量编码为整数
hitters['League'] = pd.factorize(hitters.League)[0]
hitters['Division'] = pd.factorize(hitters.Division)[0]
hitters['NewLeague'] = pd.factorize(hitters.NewLeague)[0]
hitters.head()

# 让情节出现在笔记本中
import matplotlib.pyplot as plt

#“年薪”与“工资”着色点之间的散点图
hitters.plot(kind='scatter', x='Years', y='Hits', c='Salary', colormap='jet', xlim=(0, 25), ylim=(0, 250))

#定义特征：排除职业统计（以“C”开头）和回应（薪水）
feature_cols = hitters.columns[hitters.columns.str.startswith('C') == False].drop('Salary')
feature_cols

# 定义X和y
X = hitters[feature_cols]
y = hitters.Salary

###用决策树预测工资
#使用交叉验证为决策树找到最好的** max_depth **：

#为max_depth尝试的值列表
max_depth_range = range(1, 21)

# 列表来存储每个max_depth值的平均RMSE
RMSE_scores = []

# 对每个max_depth值使用10次交叉验证
from sklearn.cross_validation import cross_val_score
for depth in max_depth_range:
    treereg = DecisionTreeRegressor(max_depth=depth, random_state=1)
    MSE_scores = cross_val_score(treereg, X, y, cv=10, scoring='neg_mean_squared_error')
    RMSE_scores.append(np.mean(np.sqrt(-MSE_scores)))

# plot max_depth（x轴）对RMSE（y轴）
plt.plot(max_depth_range, RMSE_scores)
plt.xlabel('max_depth')
plt.ylabel('RMSE (lower is better)')

# 显示最好的RMSE和相应的max_depth
sorted(zip(RMSE_scores, max_depth_range))[0]

# max_depth = 2是最好的，所以适合使用该参数的树
treereg = DecisionTreeRegressor(max_depth=2, random_state=1)
treereg.fit(X, y)

#计算功能重要性
pd.DataFrame({'feature':feature_cols, 'importance':treereg.feature_importances_}).sort_values('importance')

# ## 用随机森林预测工资

from sklearn.ensemble import RandomForestRegressor
rfreg = RandomForestRegressor()
rfreg

####调整n_estimators

#一个重要的调整参数是** n_estimators **，这是应该增长的树的数量。 它应该是一个足够大的值，这个错误似乎已经“稳定”了。

#为n_estimators尝试的值列表
estimator_range = range(10, 310, 10)

# 列表来存储每个n_estimators值的平均RMSE
RMSE_scores = []

# 对n_estimators的每个值使用5-fold交叉验证（WARNING：SLOW！）
for estimator in estimator_range:
    rfreg = RandomForestRegressor(n_estimators=estimator, random_state=1)
    MSE_scores = cross_val_score(rfreg, X, y, cv=5, scoring='neg_mean_squared_error')
    RMSE_scores.append(np.mean(np.sqrt(-MSE_scores)))

# 绘制n_estimators（x轴）与RMSE（y轴）
plt.plot(estimator_range, RMSE_scores)
plt.xlabel('n_estimators')
plt.ylabel('RMSE (lower is better)')

####调整max_features

#另一个重要的调整参数是** max_features **，这是在每个split应该考虑的功能的数量。

#为max_features尝试的值列表
feature_range = range(1, len(feature_cols)+1)

# 列表来存储每个max_features值的平均RMSE
RMSE_scores = []

# 对每个max_features值使用10次交叉验证（警告：慢！）
for feature in feature_range:
    rfreg = RandomForestRegressor(n_estimators=150, max_features=feature, random_state=1)
    MSE_scores = cross_val_score(rfreg, X, y, cv=10, scoring='neg_mean_squared_error')
    RMSE_scores.append(np.mean(np.sqrt(-MSE_scores)))

# 绘图max_features（x轴）与RMSE（y轴）
plt.plot(feature_range, RMSE_scores)
plt.xlabel('max_features')
plt.ylabel('RMSE (lower is better)')

#显示最好的RMSE和相应的max_features
sorted(zip(RMSE_scores, feature_range))[0]

####用最好的参数拟合随机森林

#max_features = 8是最好的，n_estimators = 150是足够大的
rfreg = RandomForestRegressor(n_estimators=150, max_features=8, oob_score=True, random_state=1)
rfreg.fit(X, y)

#计算功能重要性
pd.DataFrame({'feature':feature_cols, 'importance':rfreg.feature_importances_}).sort_values('importance')

# 计算出袋外R平方得分
rfreg.oob_score_

####减少X的最重要的功能
#检查X的形状
X.shape

第十课、高级scikit-learn

# 假数据
import pandas as pd
train = pd.DataFrame({'id':[0,1,2], 'length':[0.9,0.3,0.6], 'mass':[0.1,0.2,0.8], 'rings':[40,50,60]})
test = pd.DataFrame({'length':[0.59], 'mass':[0.79], 'rings':[54]})

# 训练数据
train

# 测试数据
test

# 定义X和y
feature_cols = ['length', 'mass', 'rings']
X = train[feature_cols]
y = train.id

# KN = K = 1
from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=1)
knn.fit(X, y)

# 它预测什么“应该”？
knn.predict(test)

# 让情节出现在笔记本中
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 14
plt.rcParams['figure.figsize'] = (5, 5)

# 为绘图创建一个“颜色”数组
import numpy as np
colors = np.array(['red', 'green', 'blue'])

# 训练数据的散点图，由id着色（0 =红色，1 =绿色，2 =蓝色）
plt.scatter(train.mass, train.rings, c=colors[train.id], s=50)

# 测试数据
plt.scatter(test.mass, test.rings, c='white', s=50)

# 添加标签
plt.xlabel('mass')
plt.ylabel('rings')
plt.title('How we interpret the data')

# 调整x限制
plt.scatter(train.mass, train.rings, c=colors[train.id], s=50)
plt.scatter(test.mass, test.rings, c='white', s=50)
plt.xlabel('mass')
plt.ylabel('rings')
plt.title('How KNN interprets the data')
plt.xlim(0, 30)

### StandardScaler如何解决问题？
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(X)
X_scaled = scaler.transform(X)

# 原始值
X.values

# 标准化价值
X_scaled

# 弄清楚它是如何标准化的
print scaler.mean_
print scaler.scale_

#手动标准化
(X.values - scaler.mean_) / scaler.scale_

###将StandardScaler应用于真实的数据集
url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data';;
col_names = ['label', 'color', 'proline']
wine = pd.read_csv(url, header=None, names=col_names, usecols=[0, 10, 13])

wine.head()

wine.describe()

# define X and y
feature_cols = ['color', 'proline']
X = wine[feature_cols]
y = wine.label

# 分成训练和测试集
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

# 标准化X_train
scaler.fit(X_train)
X_train_scaled = scaler.transform(X_train)

# 检查它是否正确标准化
print X_train_scaled[:, 0].mean()
print X_train_scaled[:, 0].std()
print X_train_scaled[:, 1].mean()
print X_train_scaled[:, 1].std()

#规范X_test
X_test_scaled = scaler.transform(X_test)

# 这是正确的吗？
print X_test_scaled[:, 0].mean()
print X_test_scaled[:, 0].std()
print X_test_scaled[:, 1].mean()
print X_test_scaled[:, 1].std()

# KNN对原始数据的准确性
knn = KNeighborsClassifier(n_neighbors=3)
knn.fit(X_train, y_train)
y_pred_class = knn.predict(X_test)
from sklearn import metrics
print metrics.accuracy_score(y_test, y_pred_class)

# KNN在缩放数据上的准确性
knn.fit(X_train_scaled, y_train)
y_pred_class = knn.predict(X_test_scaled)
print metrics.accuracy_score(y_test, y_pred_class)

# define X and y
feature_cols = ['color', 'proline']
X = wine[feature_cols]
y = wine.label

#对原始（未缩放的）数据进行适当的交叉验证
knn = KNeighborsClassifier(n_neighbors=3)
from sklearn.cross_validation import cross_val_score
cross_val_score(knn, X, y, cv=5, scoring='accuracy').mean()

# 为什么在缩放的数据上这个不正确的交叉验证？
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
cross_val_score(knn, X_scaled, y, cv=5, scoring='accuracy').mean()

#使用管道修复交叉验证过程
from sklearn.pipeline import make_pipeline
pipe = make_pipeline(StandardScaler(), KNeighborsClassifier(n_neighbors=3))
cross_val_score(pipe, X, y, cv=5, scoring='accuracy').mean()

#使用GridSearchCV搜索最佳的n_neighbors值
neighbors_range = range(1, 21)
param_grid = dict(kneighborsclassifier__n_neighbors=neighbors_range)
from sklearn.grid_search import GridSearchCV
grid = GridSearchCV(pipe, param_grid, cv=5, scoring='accuracy')
grid.fit(X, y)
print grid.best_score_
print grid.best_params_

第十一课、集群

#1. K-均值聚类
#2.集群评估
#3. DBSCAN集群

# 啤酒数据集
import pandas as pd
url = 'https://raw.githubusercontent.com/justmarkham/DAT8/master/data/beer.txt';;
beer = pd.read_csv(url, sep=' ')
beer

# 你将如何聚集这些啤酒？

# define X
X = beer.drop('name', axis=1)

##第一部分：K均值聚类

#有3个集群的K-means
from sklearn.cluster import KMeans
km = KMeans(n_clusters=3, random_state=1)
km.fit(X)

# 查看集群标签
km.labels_

# 保存群集标签并按群集排序
beer['cluster'] = km.labels_
beer.sort_values('cluster')

# 审查集群中心
km.cluster_centers_

# 计算每个群集的每个特征的平均值
beer.groupby('cluster').mean()

# 保存集群中心的DataFrame
centers = beer.groupby('cluster').mean()

import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 14

#为绘图创建一个“颜色”数组
import numpy as np
colors = np.array(['red', 'green', 'blue', 'yellow'])

#卡路里与酒精的散点图，由簇着色（0 =红色，1 =绿色，2 =蓝色）
plt.scatter(beer.calories, beer.alcohol, c=colors[beer.cluster], s=50)

# 由“+”标记的聚类中心
plt.scatter(centers.calories, centers.alcohol, linewidths=3, marker='+', s=300, c='black')

# 添加标签
plt.xlabel('calories')
plt.ylabel('alcohol')

# 散点图矩阵（0 =红色，1 =绿色，2 =蓝色）
pd.scatter_matrix(X, c=colors[beer.cluster], figsize=(10,10), s=100)

# ### 重复缩放数据

# 中心和缩放数据
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# K-means在缩放数据上具有3个群集
km = KMeans(n_clusters=3, random_state=1)
km.fit(X_scaled)

# 保存群集标签并按群集排序
beer['cluster'] = km.labels_
beer.sort_values('cluster')

# 审查集群中心
beer.groupby('cluster').mean()

# （0 =红色，1 =绿色，2 =蓝色）的散点图矩阵
pd.scatter_matrix(X, c=colors[beer.cluster], figsize=(10,10), s=100)

##第二部分：聚类评估
# 计算SC为K = 3
from sklearn import metrics
metrics.silhouette_score(X_scaled, km.labels_)

# 计算K = 2到K = 19的SC
k_range = range(2, 20)
scores = []
for k in k_range:
    km = KMeans(n_clusters=k, random_state=1)
    km.fit(X_scaled)
    scores.append(metrics.silhouette_score(X_scaled, km.labels_))

#绘制结果
plt.plot(k_range, scores)
plt.xlabel('Number of clusters')
plt.ylabel('Silhouette Coefficient')
plt.grid(True)

# K-means在缩放数据上具有4个群集
km = KMeans(n_clusters=4, random_state=1)
km.fit(X_scaled)
beer['cluster'] = km.labels_
beer.sort_values('cluster')

##第3部分：DBSCAN聚类

#DBSCAN，eps = 1和min_samples = 3
from sklearn.cluster import DBSCAN
db = DBSCAN(eps=1, min_samples=3)
db.fit(X_scaled)

# 查看集群标签
db.labels_

# 保存群集标签并按群集排序
beer['cluster'] = db.labels_
beer.sort_values('cluster')

# 审查集群中心
beer.groupby('cluster').mean()

# DBSCAN聚类分配的散点图矩阵（0 =红色，1 =绿色，2 =蓝色，-1 =黄色）
pd.scatter_matrix(X, c=colors[beer.cluster], figsize=(10,10), s=100)

第十二课、正则化线性模型

# 正则化
# 1.过度配合（复习）
# 2.过度拟合线性模型
# 3.线性模型的正则化
# 4.在scikit-learn中正则化回归
# 5. scikit-learn中的规范化分类
# 6.将正则化的线性模型与非正则化的线性模型进行比较

# ###加载并准备犯罪数据集

import pandas as pd
url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/communities/communities.data';;
crime = pd.read_csv(url, header=None, na_values=['?'])
crime.head()

# 检查响应变量
crime[127].describe()

# 删除分类功能
crime.drop([0, 1, 2, 3, 4], axis=1, inplace=True)

# 删除任何缺少值的行
crime.dropna(inplace=True)

# 检查形状
crime.shape

# 定义X和y
X = crime.drop(127, axis=1)
y = crime[127]

# 分成训练和测试集
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

####线性回归

#建立一个线性回归模型
from sklearn.linear_model import LinearRegression
linreg = LinearRegression()
linreg.fit(X_train, y_train)

# 检查系数
print linreg.coef_

#作出预测
y_pred = linreg.predict(X_test)

# 计算RMSE
from sklearn import metrics
import numpy as np
print np.sqrt(metrics.mean_squared_error(y_test, y_pred))

###岭回归
from sklearn.linear_model import Ridge
ridgereg = Ridge(alpha=0, normalize=True)
ridgereg.fit(X_train, y_train)
y_pred = ridgereg.predict(X_test)
print np.sqrt(metrics.mean_squared_error(y_test, y_pred))

# try alpha=0.1
ridgereg = Ridge(alpha=0.1, normalize=True)
ridgereg.fit(X_train, y_train)
y_pred = ridgereg.predict(X_test)
print np.sqrt(metrics.mean_squared_error(y_test, y_pred))

# 检查系数
print ridgereg.coef_

# 创建一个alpha值的数组
alpha_range = 10.**np.arange(-2, 3)
alpha_range

# 用RidgeCV选择最好的alpha
from sklearn.linear_model import RidgeCV
ridgeregcv = RidgeCV(alphas=alpha_range, normalize=True, scoring='mean_squared_error')
ridgeregcv.fit(X_train, y_train)
ridgeregcv.alpha_

# 预测方法使用最佳的alpha值
y_pred = ridgeregcv.predict(X_test)
print np.sqrt(metrics.mean_squared_error(y_test, y_pred))

### Lasso回归

# try alpha=0.001 and examine coefficients
from sklearn.linear_model import Lasso
lassoreg = Lasso(alpha=0.001, normalize=True)
lassoreg.fit(X_train, y_train)
print lassoreg.coef_

# 尝试α= 0.01并检查系数
lassoreg = Lasso(alpha=0.01, normalize=True)
lassoreg.fit(X_train, y_train)
print lassoreg.coef_

# 计算RMSE（对于α= 0.01）
y_pred = lassoreg.predict(X_test)
print np.sqrt(metrics.mean_squared_error(y_test, y_pred))

# 用LassoCV选择最好的alpha
from sklearn.linear_model import LassoCV
lassoregcv = LassoCV(n_alphas=100, normalize=True, random_state=1)
lassoregcv.fit(X_train, y_train)
lassoregcv.alpha_

# 检查系数
print lassoregcv.coef_

# 预测方法使用最佳的alpha值
y_pred = lassoregcv.predict(X_test)
print np.sqrt(metrics.mean_squared_error(y_test, y_pred))

##第五部分：scikit-learn中的规范化分类

# read in the dataset
url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data';;
wine = pd.read_csv(url, header=None)
wine.head()

#检查响应变量
wine[0].value_counts()

# define X and y
X = wine.drop(0, axis=1)
y = wine[0]

#分成训练和测试集
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

###逻辑回归（非正规化）
# 建立逻辑回归模型
from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression(C=1e9)
logreg.fit(X_train, y_train)

# 检查系数
print logreg.coef_

# 生成预测概率
y_pred_prob = logreg.predict_proba(X_test)
print y_pred_prob

# 计算对数损失
print metrics.log_loss(y_test, y_pred_prob)

# ###逻辑回归（正则化）
# 标准化X_train和X_test
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(X_train)
X_train_scaled = scaler.transform(X_train)
X_test_scaled = scaler.transform(X_test)

# 尝试C = 0.1与L1罚款
logreg = LogisticRegression(C=0.1, penalty='l1')
logreg.fit(X_train_scaled, y_train)
print logreg.coef_

# 生成预测概率并计算对数损失
y_pred_prob = logreg.predict_proba(X_test_scaled)
print metrics.log_loss(y_test, y_pred_prob)

# 尝试C = 0.1与二级惩罚
logreg = LogisticRegression(C=0.1, penalty='l2')
logreg.fit(X_train_scaled, y_train)
print logreg.coef_

# 生成预测概率并计算对数损失
y_pred_prob = logreg.predict_proba(X_test_scaled)
print metrics.log_loss(y_test, y_pred_prob)

#StandardScaler和LogisticRegression的管道
from sklearn.pipeline import make_pipeline
pipe = make_pipeline(StandardScaler(), LogisticRegression())

# 网格搜索C和惩罚的最佳组合
from sklearn.grid_search import GridSearchCV
C_range = 10.**np.arange(-2, 3)
penalty_options = ['l1', 'l2']
param_grid = dict(logisticregression__C=C_range, logisticregression__penalty=penalty_options)
grid = GridSearchCV(pipe, param_grid, cv=10, scoring='neg_log_loss')
grid.fit(X, y)

# 打印所有日志丢失分数
grid.grid_scores_

# 检查最好的模型
print grid.best_score_
print grid.best_params_

###第6部分：比较正规化的线性模型与未经调整的线性模型
# **正则化线性模型的优点：**
#  - 更好的性能
# L1正则化执行自动特征选择
#  - 用于高维问题（p> n）
# **正则化线性模型的缺点：**
#  - 调整是必需的
#  - 推荐功能缩放
#  - 较少可解释（由于功能缩放）

第五课、 比较多项式和高斯朴素贝叶斯

# read the data
import pandas as pd
url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/pima-indians-diabetes/pima-indians-diabetes.data';;
col_names = ['pregnant', 'glucose', 'bp', 'skin', 'insulin', 'bmi', 'pedigree', 'age', 'label']
pima = pd.read_csv(url, header=None, names=col_names)

#注意所有功能都是连续的
pima.head()

# 创建X和Y.
X = pima.drop('label', axis=1)
y = pima.label

# 分成训练和测试集
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

# 导入多项式和高斯朴素贝叶斯
from sklearn.naive_bayes import MultinomialNB, GaussianNB
from sklearn import metrics

# 多项式朴素贝叶斯检验的准确性
mnb = MultinomialNB()
mnb.fit(X_train, y_train)
y_pred_class = mnb.predict(X_test)
print metrics.accuracy_score(y_test, y_pred_class)

# 高斯朴素贝叶斯检验的准确性
gnb = GaussianNB()
gnb.fit(X_train, y_train)
y_pred_class = gnb.predict(X_test)
print metrics.accuracy_score(y_test, y_pred_class)

#**结论：**将朴素贝叶斯分类应用于具有**连续特征**的数据集时，最好使用高斯朴素贝叶斯比多项朴素贝叶斯。
# 后者适用于包含**离散特征**（例如，字数）的数据集。

第六课、自然语言处理（NLP）

#NLP需要理解**语言**和**世界**。
###第一部分：阅读Yelp评论
# - “语料库”=文件的集合
# - “corpora”= corpu的复数形式s

import pandas as pd
import numpy as np
import scipy as sp
from sklearn.cross_validation import train_test_split
from sklearn.feature_extraction.text import CountVectorizer, TfidfVectorizer
from sklearn.naive_bayes import MultinomialNB
from sklearn.linear_model import LogisticRegression
from sklearn import metrics
from textblob import TextBlob, Word
from nltk.stem.snowball import SnowballStemmer

# 将yelp.csv读入DataFrame
url = 'https://raw.githubusercontent.com/justmarkham/DAT8/master/data/yelp.csv';;
yelp = pd.read_csv(url)

# 创建一个只包含五星评论和一星评论的新数据框
yelp_best_worst = yelp[(yelp.stars == 5) | (yelp.stars == 1)]

# 定义X和y
X = yelp_best_worst.text
y = yelp_best_worst.stars

# 将新的DataFrame分成训练和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

###第2部分：令牌化

# - **什么：**将文本分成单词，如句子或单词
# - **为什么：**给以前非结构化文本的结构
# - **备注：**英文文本比较容易，一些语言不容易

#使用CountVectorizer从X_train和X_test创建文档项矩阵
vect = CountVectorizer()
X_train_dtm = vect.fit_transform(X_train)
X_test_dtm = vect.transform(X_test)

# 行是文档，列是术语（又名“令牌”或“功能”）
X_train_dtm.shape

# 最后50个功能
print vect.get_feature_names()[-50:]

#显示矢量选项
vect

# - **小写：**布尔值，默认为true
# - 在标记之前将所有字符转换为小写。

# 不要转换为小写
vect = CountVectorizer(lowercase=False)
X_train_dtm = vect.fit_transform(X_train)
X_train_dtm.shape

# - ** ngram_range：**元组（min_n，max_n）
# - 要提取的不同n-gram的n值范围的下限和上限。 将使用n的所有值，使得min_n <= n <= max_n。
#包括1克和2克

vect = CountVectorizer(ngram_range=(1, 2))
X_train_dtm = vect.fit_transform(X_train)
X_train_dtm.shape

# 最后50个功能
print vect.get_feature_names()[-50:]

#**预测星级评分：**

#使用CountVectorizer的默认选项
vect = CountVectorizer()

# 创建文档术语矩阵
X_train_dtm = vect.fit_transform(X_train)
X_test_dtm = vect.transform(X_test)

# 使用朴素贝叶斯来预测星级
nb = MultinomialNB()
nb.fit(X_train_dtm, y_train)
y_pred_class = nb.predict(X_test_dtm)

# 计算准确度
print metrics.accuracy_score(y_test, y_pred_class)

# 计算零准确度
y_test_binary = np.where(y_test == 5, 1, 0)
max(y_test_binary.mean(), 1 - y_test_binary.mean())

# 定义一个接受向量的函数并计算精度
def tokenize_test(vect):
    X_train_dtm = vect.fit_transform(X_train)
    print 'Features: ', X_train_dtm.shape[1]
    X_test_dtm = vect.transform(X_test)
    nb = MultinomialNB()
    nb.fit(X_train_dtm, y_train)
    y_pred_class = nb.predict(X_test_dtm)
    print 'Accuracy: ', metrics.accuracy_score(y_test, y_pred_class)

# 包括1克和2克
vect = CountVectorizer(ngram_range=(1, 2))
tokenize_test(vect)

##第3部分：停用词删除
# - **什么：**删除可能出现在任何文本中的常见词汇
# - **为什么：**他们没有告诉你很多关于你的文字
#show vectorizer选项

vect

# - ** stop_words：** string {'english'}，列表或None（默认）
# - 如果是英文，则使用英文的内置停用词表。
# - 如果一个列表，假设该列表包含停用词，所有这些都将从结果标记中删除。
# - 如果没有，则不会使用停用词。 可以将max_df设置为[0.7,1.0]范围内的一个值，以便根据术语内部语料文档的频率自动检测并过滤停用词。
# 删除英文停用词
vect = CountVectorizer(stop_words='english')
tokenize_test(vect)

# 停止词的集合
print vect.get_stop_words()

###第4部分：其他CountVectorizer选项

# - ** max_features：** int或None，default = None
# - 如果不是None，建立一个词汇表，只考虑整个语料库中词频排序的顶级max_features。
#删除英文停用词，只保留100个功能
vect = CountVectorizer(stop_words='english', max_features=100)
tokenize_test(vect)

#全部100个功能
print vect.get_feature_names()

# 包括1克和2克，并限制功能的数量
vect = CountVectorizer(ngram_range=(1, 2), max_features=100000)
tokenize_test(vect)

#  - ** min_df：**浮动范围[0.0，1.0]或int，默认= 1
#  - 建立词汇表时忽略文档频率严格低于给定阈值的词汇。 这个值在文献中也被称为切断。 如果是浮点数，则该参数表示文档的一个比例，即整数绝对计数。
#
# 包括1克和2克，只包括出现至少2次的术语
vect = CountVectorizer(ngram_range=(1, 2), min_df=2)
tokenize_test(vect)

##第5部分：TextBlob简介

#TextBlob：“简体文字处理”
#打印第一个评论
print yelp_best_worst.text[0]

# 将其另存为TextBlob对象
# curl https://raw.github.com/sloria/TextBlob/master/download_corpora.py | python
#python -m textblob.download_corpora 安装语料库
review = TextBlob(yelp_best_worst.text[0])

#列出单词
review.words

# 列出句子
review.sentences

# 一些字符串方法可用
review.lower()

###第6部分：词干和词形化
# ＃**词干：**
#  - **什么：**减少一个词的基础/干/根的形式
#  - **为什么：**以相同的方式处理相关的词常常是有意义的
#  - **备注：**
#  - 使用“简单”和快速的基于规则的方法
#  - 词干通常不会显示给用户（用于分析/索引）
#  - 一些搜索引擎将同词干的单词视为同义词
#
# ＃初始化stemmer
stemmer = SnowballStemmer('english')

# 干每个字
print [stemmer.stem(word) for word in review.words]

# **词形化**
# - **什么：**推导一个单词的规范形式（“引理”）
# - **为什么：**可以比干扰更好
# - **注释：**使用基于字典的方法（比词干慢）
# ＃假设每个单词都是一个名词
print [word.lemmatize() for word in review.words]

# 假定每个单词都是动词
print [word.lemmatize(pos='v') for word in review.words]

# 定义一个接受文本并返回引理列表的函数
def split_into_lemmas(text):
    text = unicode(text, 'utf-8').lower()
    words = TextBlob(text).words
    return [word.lemmatize() for word in words]

#使用split_into_lemmas作为特征提取功能（WARNING：SLOW！）
vect = CountVectorizer(analyzer=split_into_lemmas)
tokenize_test(vect)

# 最后50个功能
print vect.get_feature_names()[-50:]

###第7部分：术语频率逆文档频率（TF-IDF）

# - **什么：**计算一个单词出现在文档中的“相对频率”，而不是所有文档的频率
# - **为什么：**比“词频”更有用于识别每个文档中的“重要”单词（该文档中频率高，其他文档中频率低）
# - **备注：**用于搜索引擎评分，文字摘要，文档聚类

# 示例文档
simple_train = ['call you tonight', 'Call me a cab', 'please call me... PLEASE!']

# 术语频率
vect = CountVectorizer()
tf = pd.DataFrame(vect.fit_transform(simple_train).toarray(), columns=vect.get_feature_names())
tf

#文档频率
vect = CountVectorizer(binary=True)
df = vect.fit_transform(simple_train).toarray().sum(axis=0)
pd.DataFrame(df.reshape(1, 6), columns=vect.get_feature_names())

# 术语频率反向文档频率（简单版本）
tf / df

# TfidfVectorizer
vect = TfidfVectorizer()
pd.DataFrame(vect.fit_transform(simple_train).toarray(), columns=vect.get_feature_names())

#**更多细节：** [TF-IDF是关于什么的]（http://planspace.org/20150524-tfidf_is_about_what_matters/）

###第8部分：使用TF-IDF总结Yelp评论
#Reddit的autotldr使用基于TF-IDF的[SMMRY]（http://smmry.com/about）算法！

# 使用TF-IDF创建一个文档项矩阵
vect = TfidfVectorizer(stop_words='english')
dtm = vect.fit_transform(yelp.text)
features = vect.get_feature_names()
dtm.shape

def summarize():
    #选择一个至少300个字符的随机评论
    review_length = 0
    while review_length < 300:
        review_id = np.random.randint(0, len(yelp))
        review_text = unicode(yelp.text[review_id], 'utf-8')
        review_length = len(review_text)

    # 创建一个单词字典和他们的TF-IDF分数
    word_scores = {}
    for word in TextBlob(review_text).words:
        word = word.lower()
        if word in features:
            word_scores[word] = dtm[review_id, features.index(word)]

    # 打印前5名TF-IDF分数的单词
    print 'TOP SCORING WORDS:'
    top_scores = sorted(word_scores.items(), key=lambda x: x[1], reverse=True)[:5]
    for word, score in top_scores:
        print word

    # 随机打印5个单词
    print '\n' + 'RANDOM WORDS:'
    random_words = np.random.choice(word_scores.keys(), size=5, replace=False)
    for word in random_words:
        print word

    # 打印评论
    print '\n' + review_text

summarize()

# ## 第九部分：情感分析

print review

# 极性范围从-1（最负）到1（最正）
review.sentiment.polarity

# 了解应用方法
yelp['length'] = yelp.text.apply(len)
yelp.head(1)

# 定义一个接受文本并返回极性的函数
def detect_sentiment(text):
    return TextBlob(text.decode('utf-8')).sentiment.polarity

# 为情绪创建一个新的DataFrame列（WARNING：SLOW！）
yelp['sentiment'] = yelp.text.apply(detect_sentiment)

# 按星星分组的情节盒情节
yelp.boxplot(column='sentiment', by='stars')

# 评论与最积极的情绪
yelp[yelp.sentiment == 1].text.head()

# 评论与最消极的情绪
yelp[yelp.sentiment == -1].text.head()

# 扩大列显示
pd.set_option('max_colwidth', 500)

# 5星评论中的负面情绪
yelp[(yelp.stars == 5) & (yelp.sentiment < -0.3)].head(1)

# 积极的情绪在一星评论
yelp[(yelp.stars == 1) & (yelp.sentiment > 0.5)].head(1)

# 重置列显示宽度
pd.reset_option('max_colwidth')

###奖金：将特征添加到文档术语矩阵

#创建一个仅包含5星评论和1星评论的DataFrame
yelp_best_worst = yelp[(yelp.stars == 5) | (yelp.stars == 1)]

# 定义X和y
feature_cols = ['text', 'sentiment', 'cool', 'useful', 'funny']
X = yelp_best_worst[feature_cols]
y = yelp_best_worst.stars

# 分成训练和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

# 仅使用带有文本列的CountVectorizer
vect = CountVectorizer()
X_train_dtm = vect.fit_transform(X_train.text)
X_test_dtm = vect.transform(X_test.text)
print X_train_dtm.shape
print X_test_dtm.shape

# 其他四个特征栏的形状
X_train.drop('text', axis=1).shape

# 将其他特征列转换为浮点并转换为稀疏矩阵
extra = sp.sparse.csr_matrix(X_train.drop('text', axis=1).astype(float))
extra.shape

# 组合稀疏矩阵
X_train_dtm_extra = sp.sparse.hstack((X_train_dtm, extra))
X_train_dtm_extra.shape

# 重复测试设置
extra = sp.sparse.csr_matrix(X_test.drop('text', axis=1).astype(float))
X_test_dtm_extra = sp.sparse.hstack((X_test_dtm, extra))
X_test_dtm_extra.shape

# 仅使用文本列进行逻辑回归
logreg = LogisticRegression(C=1e9)
logreg.fit(X_train_dtm, y_train)
y_pred_class = logreg.predict(X_test_dtm)
print metrics.accuracy_score(y_test, y_pred_class)

# 使用具有所有特征的逻辑回归
logreg = LogisticRegression(C=1e9)
logreg.fit(X_train_dtm_extra, y_train)
y_pred_class = logreg.predict(X_test_dtm_extra)
print metrics.accuracy_score(y_test, y_pred_class)

# ##奖金：有趣的TextBlob功能

# 拼写更正
TextBlob('15 minuets late').correct()

# 拼写检查
Word('parot').spellcheck()

# 定义
Word('bank').define('v')

# 语言识别
a = TextBlob('Hola amigos').detect_language()
a

###结论
#NLP是一个巨大的领域
# - 了解基础知识将扩大您可以使用的数据类型
# - 简单的技术有很长的路要走
# - 尽可能使用scikit-learn来进行NLP

第七课、用自行车数据进行练习

import pandas as pd
import numpy as np
from sklearn.cross_validation import cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.tree import DecisionTreeRegressor, export_graphviz

# 读取数据并设置“datetime”作为索引
url = 'https://raw.githubusercontent.com/justmarkham/DAT8/master/data/bikeshare.csv';;
bikes = pd.read_csv(url, index_col='datetime', parse_dates=True)

# “count”是一种方法，所以最好重命名该列
bikes.rename(columns={'count':'total'}, inplace=True)

# 创建“小时”作为自己的功能
bikes['hour'] = bikes.index.hour

bikes.head()

bikes.tail()

# ＃＃ 任务1
#运行这两个`groupby`语句，并找出他们告诉你的数据。
#“每日工作日”的平均租金
bikes.groupby('workingday').total.mean()

# 每个“小时”值的平均租金
bikes.groupby('hour').total.mean()

###任务2
#运行这个绘图代码，并确保你了解输出。 然后，把这个情节分成两个独立的情节，以“工作日”为条件。 （换句话说，一个小区应显示“workingday = 0”的小时趋势，另一个应显示“workingday = 1”的小时趋势。）

#平均租金为“小时”的每个值
bikes.groupby('hour').total.mean().plot()

# “工作日= 0”的小时租金趋势
bikes[bikes.workingday==0].groupby('hour').total.mean().plot()

# “工作日= 1”每小时租金趋势
bikes[bikes.workingday==1].groupby('hour').total.mean().plot()

# 结合这两块地块
bikes.groupby(['hour', 'workingday']).total.mean().unstack().plot()

###任务3
#使用“total”作为响应，“hour”和“workingday”作为唯一的特征，将线性回归模型拟合到整个数据集中。 然后，打印系数并解释它们。 在这种情况下线性回归的局限性是什么？
#创建X和Y
feature_cols = ['hour', 'workingday']
X = bikes[feature_cols]
y = bikes.total

# 拟合线性回归模型并打印系数
linreg = LinearRegression()
linreg.fit(X, y)
linreg.coef_

###任务4

#使用10倍交叉验证来计算线性回归模型的RMSE。

#保存10个由cross_val_score输出的MSE分数
scores = cross_val_score(linreg, X, y, cv=10, scoring='neg_mean_squared_error')

# 将MSE转换为RMSE，然后计算10个RMSE分数的均值
np.mean(np.sqrt(-scores))

###任务5
#使用10倍交叉验证来评估具有相同功能的决策树模型（适合您选择的任何“max_depth”）。
#用“max_depth = 7”评估决策树模型

treereg = DecisionTreeRegressor(max_depth=7, random_state=1)
scores = cross_val_score(treereg, X, y, cv=10, scoring='neg_mean_squared_error')
np.mean(np.sqrt(-scores))

###任务6
#使用“max_depth = 3”将决策树模型拟合到整个数据集中，并使用Graphviz创建树图。 然后，找出每个叶代表什么。 决策树学到一个线性回归模型不能学习的是什么？
#用“max_depth = 3”来拟合决策树模型
treereg = DecisionTreeRegressor(max_depth=3, random_state=1)
treereg.fit(X, y)

# 创建一个Graphviz文件
export_graphviz(treereg, out_file='tree_bikeshare.dot', feature_names=feature_cols)

#在命令行运行这个转换为PNG：
#dot -Tpng tree_bikeshare.dot -o tree_bikeshare.pn
#！[自行车数据树]（images / tree_bikeshare.png）

第一课、高级模型评估

# 第一部分：处理缺失值
import pandas as pd
url = 'https://raw.githubusercontent.com/justmarkham/DAT8/master/data/titanic.csv';;
titanic = pd.read_csv(url, index_col='PassengerId')
titanic.shape

# 检查缺少的值
titanic.isnull().sum()

#一个可能的策略是**丢失缺失值**：

#删除任何缺少值的行
titanic.dropna().shape

# 删除年龄缺失的行
titanic[titanic.Age.notnull()].shape

#有时候更好的策略是**推测缺失值**：

# 平均年龄
titanic.Age.mean()

#中年的
titanic.Age.median()

# 最常见的年龄
titanic.Age.mode()

# 填写年龄中位数的缺失值
titanic.Age.fillna(titanic.Age.median(), inplace=True)

# ##第2部分：处理分类特征（复审）

titanic.head(10)

# 编码Sex_Female功能
titanic['Sex_Female'] = titanic.Sex.map({'male':0, 'female':1})

# 创建一个虚拟变量的数据框
embarked_dummies = pd.get_dummies(titanic.Embarked, prefix='Embarked')
embarked_dummies.drop(embarked_dummies.columns[0], axis=1, inplace=True)

# 连接原始DataFrame和虚拟DataFrame
titanic = pd.concat([titanic, embarked_dummies], axis=1)

titanic.head(1)

# 定义X和y
feature_cols = ['Pclass', 'Parch', 'Age', 'Sex_Female', 'Embarked_Q', 'Embarked_S']
X = titanic[feature_cols]
y = titanic.Survived

# train/test 分割
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

# 训练逻辑回归模型
from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression(C=1e9)
logreg.fit(X_train, y_train)

# 对测试集进行预测
y_pred_class = logreg.predict(X_test)

# 计算测试精度
from sklearn import metrics
print metrics.accuracy_score(y_test, y_pred_class)

# ## 第3部分：ROC曲线和AUC

# 预测生存概率
y_pred_prob = logreg.predict_proba(X_test)[:, 1]

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (8, 6)
plt.rcParams['font.size'] = 14

# 绘制ROC曲线
# ROC曲线指受试者工作特征曲线 / 接收器操作特性曲线(receiver operating characteristic curve), \
# 是反映敏感性和特异性连续变量的综合指标,是用构图法揭示敏感性和特异性的相互关系，
# 它通过将连续变量设定出多个不同的临界值，从而计算出一系列敏感性和特异性，再以敏感性为纵坐标、（1-特异性）为横坐标绘制成曲线，
# 曲线下面积越大，诊断准确性越高。在ROC曲线上，最靠近坐标图左上方的点为敏感性和特异性均较高的临界值。

fpr, tpr, thresholds = metrics.roc_curve(y_test, y_pred_prob)
plt.plot(fpr, tpr)
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.0])
plt.xlabel('False Positive Rate (1 - Specificity)')
plt.ylabel('True Positive Rate (Sensitivity)')

# 计算AUC
# AUC（Area Under Curve）被定义为ROC曲线下的面积，显然这个面积的数值不会大于1。
# 又由于ROC曲线一般都处于y=x这条直线的上方，所以AUC的取值范围在0.5和1之间。
# 使用AUC值作为评价标准是因为很多时候ROC曲线并不能清晰的说明哪个分类器的效果更好，而作为一个数值，对应AUC更大的分类器效果更好。

print metrics.roc_auc_score(y_test, y_pred_prob)

#除了允许您计算AUC之外，查看ROC曲线可以帮助您选择一个阈值，**以平衡灵敏度和特异性**的方式对特定上下文有意义。
#根据实际响应值分组的预测概率直方图
df = pd.DataFrame({'probability':y_pred_prob, 'actual':y_test})
df.hist(column='probability', by='actual', sharex=True, sharey=True)

#如果您在绘制ROC曲线或计算AUC时使用了** y_pred_class **而不是** y_pred_prob **，会发生什么？
# 使用y_pred_class的ROC曲线 - 错误！
fpr, tpr, thresholds = metrics.roc_curve(y_test, y_pred_class)
plt.plot(fpr, tpr)

# AUC使用y_pred_class - 错误！
print metrics.roc_auc_score(y_test, y_pred_class)

#如果使用** y_pred_class **，它将把0和100解释为预测概率0％和100％。
## ##奖金：ROC曲线只对预测概率的排名顺序敏
#打印前10个预测概率
y_pred_prob[:10]

# 取预测概率的平方根（使它们全部变大）
import numpy as np
y_pred_prob_new = np.sqrt(y_pred_prob)

# 打印修改后的预测概率
y_pred_prob_new[:10]

# 预测概率的直方图已经改变
df = pd.DataFrame({'probability':y_pred_prob_new, 'actual':y_test})
df.hist(column='probability', by='actual', sharex=True, sharey=True)

# ROC曲线没有改变
fpr, tpr, thresholds = metrics.roc_curve(y_test, y_pred_prob_new)
plt.plot(fpr, tpr)

#AUC没有改变
print metrics.roc_auc_score(y_test, y_pred_prob_new)

##第四部分：交叉验证

#计算交叉验证的AUC
from sklearn.cross_validation import cross_val_score
cross_val_score(logreg, X, y, cv=10, scoring='roc_auc').mean()

# 添加票价模型
feature_cols = ['Pclass', 'Parch', 'Age', 'Sex_Female', 'Embarked_Q', 'Embarked_S', 'Fare']
X = titanic[feature_cols]

# 重新计算AUC
cross_val_score(logreg, X, y, cv=10, scoring='roc_auc').mean()

第二课、交叉验证参数调整，模型选择和特征选择

##目标：评估模型评估程序、需要一种在机器学习模型之间进行选择的方法
# 估计**样本外的数据的可能性能

from sklearn.datasets import load_iris
from sklearn.cross_validation import train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn import metrics

# 读入iris数据
iris = load_iris()

# 创建X（特征）和Y（响应）
X = iris.data
y = iris.target

# 使用train / test拆分不同的random_state值

# train_test_split是交叉验证中常用的函数，功能是从样本中随机的按比例选取train data和test data，形式为：

# 参数代表含义:
# train_data：所要划分的样本特征集
# train_target：所要划分的样本结果
# test_size：样本占比，如果是整数的话就是样本的数量
# random_state：是随机数的种子。

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=4)

#用K = 5检查KNN的分类精度
knn = KNeighborsClassifier(n_neighbors=5)
knn.fit(X_train, y_train)
y_pred = knn.predict(X_test)
print metrics.accuracy_score(y_test, y_pred)

# 模拟将25个观察数据集分成5个等级
from sklearn.cross_validation import KFold
kf = KFold(25, n_folds=5, shuffle=False)

# 打印每个训练和测试集的内容
print '{} {:^61} {}'.format('Iteration', 'Training set observations', 'Testing set observations')
for iteration, data in enumerate(kf, start=1):
    print '{:^9} {} {:^25}'.format(iteration, data[0], data[1])

#  - 数据集包含** 25个观察**（编号0到24）
#  - 5倍交叉验证，因此它运行** 5次迭代**
#  - 对于每一次迭代，每个观察值都在训练集或测试集**中，但不是两者**
#  - 每个观察都在测试集**中，只有一次**

###比较交叉验证到训练/测试分割
# **交叉验证的优点：**
#  - 更准确地估计样本外的准确性
#  - 更有效地使用数据（每次观察都用于培训和测试）

# **列车/测试分组的优点：**
#  - 比K-fold交叉验证运行速度快K倍
#  - 更简单地检查测试过程的详细结果

##交叉验证示例：参数调整

#**目标：**为iris数据集上的KNN选择最佳调整参数（又名“超参数”）
from sklearn.cross_validation import cross_val_score

# 对于KNN（n_neighbors参数），K = 5的10倍交叉验证
knn = KNeighborsClassifier(n_neighbors=5)
scores = cross_val_score(knn, X, y, cv=10, scoring='accuracy')
print scores

# 使用平均精度作为样本外精度的估计
print scores.mean()

#为KNN搜索K的最优值
k_range = range(1, 31)
k_scores = []
for k in k_range:
    knn = KNeighborsClassifier(n_neighbors=k)
    scores = cross_val_score(knn, X, y, cv=10, scoring='accuracy')
    k_scores.append(scores.mean())
print k_scores

import matplotlib.pyplot as plt

# 绘制KNN（x轴）的K值与交叉验证的精度（y轴）
plt.plot(k_range, k_scores)
plt.xlabel('Value of K for KNN')
plt.ylabel('Cross-Validated Accuracy')

###交叉验证示例：模型选择

#**目标：**比较最佳KNN模型和iris数据集上的逻辑回归
#用最好的KNN模型进行10次交叉验证

knn = KNeighborsClassifier(n_neighbors=20)
print cross_val_score(knn, X, y, cv=10, scoring='accuracy').mean()

# 逻辑回归10倍交叉验证
from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression()
print cross_val_score(logreg, X, y, cv=10, scoring='accuracy').mean()

# ##交叉验证示例：功能选择

# ＃**目标**：选择报纸特征是否应包含在广告数据集的线性回归模型中
import pandas as pd
import numpy as np
from sklearn.linear_model import LinearRegression

# 读入广告数据集
data = pd.read_csv('http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv';;, index_col=0)

# 创建三个功能名称的Python列表
feature_cols = ['TV', 'radio', 'newspaper']   #注意key值还区分大小写

# 使用列表来选择DataFrame（X）的一个子集
X = data[feature_cols]

# 选择Sales列作为回应（y）
y = data.sales

# 所有三个功能的10倍交叉验证
lm = LinearRegression()
scores = cross_val_score(lm, X, y, cv=10, scoring='neg_mean_squared_error')
print scores

# 修复MSE分数的标志
mse_scores = -scores
print mse_scores

#从MSE转换到RMSE
rmse_scores = np.sqrt(mse_scores)
print rmse_scores

# 计算平均RMSE
print rmse_scores.mean()

# 具有两种功能的10倍交叉验证（不包括报纸）
feature_cols = ['TV', 'radio']
X = data[feature_cols]
print np.sqrt(-cross_val_score(lm, X, y, cv=10, scoring='neg_mean_squared_error')).mean()

##改进交叉验证

# **重复交叉验证**
#  - 多次重复交叉验证（使用**不同的随机分割**数据）并对结果进行平均
#  - 通过**减少与交叉验证的单个试验相关联的方差**，更可靠地估计样本外的表现
# **创建一个保留集**
# **在开始模型构建过程之前，“拿出”一部分数据**
#  - 使用剩余数据的交叉验证找到最佳模型，并使用保留集**对其进行测试**
#  - 样本外表现更可靠的估计，因为持有集**是真正的样本外**
# **交叉验证迭代中的特征工程和选择**
#  - 通常，在交叉验证之前，要素工程和选择发生**
#  - 而是在每次交叉验证迭代中执行所有特征工程和选择**
#  - 更可靠的样本外表现评估，因为它更好地模拟**模型在样本外的应用

第三课、 初识朴素贝叶斯分类

*P( 类别 | 特征) = P ( 特征 | 类别 ) P( 类别) / P(特征)**

import pandas as pd
import numpy as np

# 将iris数据读入DataFrame
url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data';;
col_names = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'species']  #萼片（sepal）、花瓣（petal） species（种类）
iris = pd.read_csv(url, header=None, names=col_names)
iris.head()

# 将天花板函数应用于数字列
iris.loc[:, 'sepal_length':'petal_width'] = iris.loc[:, 'sepal_length':'petal_width'].apply(np.ceil)
iris.head()

# ## 决定如何进行预测
#假设我有一个**样本外的光圈**，具有以下测量值：** 7，3，5，2 **。 我如何预测这个物种？

#显示所有具有以下特征的观察值：7,3,5,2
iris[(iris.sepal_length==7) & (iris.sepal_width==3) & (iris.petal_length==5) & (iris.petal_width==2)]

# 为这些观测计数物种
iris[(iris.sepal_length==7) & (iris.sepal_width==3) & (iris.petal_length==5) & (iris.petal_width==2)].species.value_counts()

# 计数tr所有观察物种fo
iris.species.value_counts()

第四课、朴素贝叶斯分类处理文本文件

from sklearn.feature_extraction.text import CountVectorizer

# 从一个简单的例子开始
simple_train = ['call you tonight', 'Call me a cab', 'please call me... PLEASE!']

# 学习训练数据的“词汇”
vect = CountVectorizer()
vect.fit(simple_train)
vect.get_feature_names()

# 将训练数据转换成“文档 - 术语矩阵”
simple_train_dtm = vect.transform(simple_train)
simple_train_dtm

#打印稀疏矩阵
print simple_train_dtm

# 将稀疏矩阵转换为稠密矩阵
simple_train_dtm.toarray()

# 一起检查词汇和文档项目矩阵
import pandas as pd
pd.DataFrame(simple_train_dtm.toarray(), columns=vect.get_feature_names())

#> - 每个单独的标记出现频率（标准化或不标准）被视为一个**功能**。
#> - 给定文档的所有标记频率的矢量被认为是一个多元**样本**。
#因此，**文档**可以由一个矩阵表示，每个文档有一行**和每个记号**有一列（例如单词）出现在语料库中。

#>我们把**矢量化称为将文本文档集合转化为数字特征向量的一般过程。这个具体的策略（标记化，计数和标准化）被称为**袋子或**袋子的代表。
# 文档通过单词出现来描述，而完全忽略文档中单词的相对位置信息。

#将测试数据转换成文档术语矩阵（使用现有的词汇）
simple_test = ["please don't call me"]
simple_test_dtm = vect.transform(simple_test)
simple_test_dtm.toarray()

# 一起检查词汇和文档项目矩阵
pd.DataFrame(simple_test_dtm.toarray(), columns=vect.get_feature_names())

#**总结：**
# - `vect.fit（train）`学习训练数据的词汇
# - `vect.transform（train）`使用拟合词汇从训练数据建立一个文档项矩阵
# - `vect.transform（test）`使用拟合的词汇从测试数据中建立一个文档项矩阵（并忽略它以前没见过的记号）
# ## Part 2: 读取SMS数据

# 读取制表符分隔的文件
url = 'https://raw.githubusercontent.com/justmarkham/DAT8/master/data/sms.tsv';;
col_names = ['label', 'message']
sms = pd.read_table(url, sep='\t', header=None, names=col_names)
print sms.shape

sms.head(20)

sms.label.value_counts()

#将标签转换为数字变量
sms['label'] = sms.label.map({'ham':0, 'spam':1})

# 定义X和y
X = sms.message
y = sms.label

#分成训练和测试集
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)
print X_train.shape
print X_test.shape

# ## Part 3: 向量化SMS数据

# 实例化矢量化器
vect = CountVectorizer()

# 学习训练数据词汇，然后创建文档术语矩阵
vect.fit(X_train)
X_train_dtm = vect.transform(X_train)
X_train_dtm

# 替代方法：将适合和变换合并为一个步骤
X_train_dtm = vect.fit_transform(X_train)
X_train_dtm

# 将测试数据（使用拟合的词汇）转换成文档术语矩阵
X_test_dtm = vect.transform(X_test)
X_test_dtm

# ## Part 4: 检查令牌和他们的计数

# 存储令牌名称
X_train_tokens = vect.get_feature_names()

# 前50个令牌
print X_train_tokens[:50]

# 最后50个令牌
print X_train_tokens[-50:]

# 将X_train_dtm视为一个稠密矩阵
X_train_dtm.toarray()

# 计算每个令牌在X_train_dtm中的所有消息中出现多少次
import numpy as np
X_train_counts = np.sum(X_train_dtm.toarray(), axis=0)
X_train_counts

X_train_counts.shape

# 用它们的计数创建一个令牌的DataFrame
pd.DataFrame({'token':X_train_tokens, 'count':X_train_counts}).sort_values('count')

# ## 计算每个令牌的“spamminess”

#为火腿和垃圾邮件创建单独的数据框
sms_ham = sms[sms.label==0]
sms_spam = sms[sms.label==1]

#学习所有消息的词汇并保存
vect.fit(sms.message)
all_tokens = vect.get_feature_names()

# 为火腿和垃圾邮件创建文档术语矩阵
ham_dtm = vect.transform(sms_ham.message)
spam_dtm = vect.transform(sms_spam.message)

# 计算EACH令牌在所有的火腿消息中出现的次数
ham_counts = np.sum(ham_dtm.toarray(), axis=0)

# 统计每个垃圾邮件中每个令牌的显示次数
spam_counts = np.sum(spam_dtm.toarray(), axis=0)

# 使用单独的火腿和垃圾邮件计数创建一个令牌数据框
token_counts = pd.DataFrame({'token':all_tokens, 'ham':ham_counts, 'spam':spam_counts})

# 添加一个火腿和垃圾邮件计数以避免被零除（在下面的步骤中）
token_counts['ham'] = token_counts.ham + 1
token_counts['spam'] = token_counts.spam + 1

#计算每个令牌的垃圾邮件到火腿的比例
token_counts['spam_ratio'] = token_counts.spam / token_counts.ham
token_counts.sort_values('spam_ratio')

# 第5部分：建立一个朴素贝叶斯模型
from sklearn.naive_bayes import MultinomialNB
nb = MultinomialNB()
nb.fit(X_train_dtm, y_train)

# 对X_test_dtm进行类别预测
y_pred_class = nb.predict(X_test_dtm)

# 计算类别预测的准确性
from sklearn import metrics
print metrics.accuracy_score(y_test, y_pred_class)

# 混乱矩阵
print metrics.confusion_matrix(y_test, y_pred_class)

# 预测（校准不良）概率
y_pred_prob = nb.predict_proba(X_test_dtm)[:, 1]
y_pred_prob

# 计算AUC
print metrics.roc_auc_score(y_test, y_pred_prob)

#打印错误信息的消息文本
X_test[y_test < y_pred_class]

# 打印消息文本的错误否定
X_test[y_test > y_pred_class]

# 你有什么注意到这个漏报？
X_test[3132]

# ## Part 6: 比较朴素贝叶斯和逻辑回归

# 进口/实例化/ FIT
from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression(C=1e9)
logreg.fit(X_train_dtm, y_train)

# 类别预测和预测概率
y_pred_class = logreg.predict(X_test_dtm)
y_pred_prob = logreg.predict_proba(X_test_dtm)[:, 1]

# 计算准确度和AUC
print metrics.accuracy_score(y_test, y_pred_class)
print metrics.roc_auc_score(y_test, y_pred_prob)

第一课： 线性回归分析bikes数据

# 读取数据并将日期时间设置为索引
import pandas as pd
url = 'https://raw.githubusercontent.com/justmarkham/DAT8/master/data/bikeshare.csv';;
bikes = pd.read_csv(url, index_col='datetime', parse_dates=True)
bikes.head()

# **问题：**
#  - 每个观察代表什么？
#  - 什么是响应变量（由Kaggle定义）？
#  - 有多少功能？
# ＃“count”是一种方法，所以最好把这个列命名为别的

bikes.rename(columns={'count':'total'}, inplace=True)

# ## 可视化数据

import seaborn as sns
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (8, 6)
plt.rcParams['font.size'] = 14

# 大熊猫散点图
bikes.plot(kind='scatter', x='temp', y='total', alpha=0.2)

# Seaborn散点图与回归线
sns.lmplot(x='temp', y='total', data=bikes, aspect=1.5, scatter_kws={'alpha':0.2})

##线性回归的形式：$y = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... + \beta_nx_n$

# - $ y $是回应
# - $ \ beta_0 $是拦截
# - $ \ beta_1 $是$ x_1 $的系数（第一项功能）
# - $ \ beta_n $是$ x_n $的系数（第n个特征）
#$ $ beta $值被称为**模型系数**：

# - 在模型拟合过程中使用**最小平方标准**估计（或“学习”）这些值。
# - 具体来说，我们找到最小化**平方残差**（或“平方误差之和”）的线（数学）。
# - 一旦我们学习了这些系数，我们可以使用模型来预测响应

#创建X和Y.
feature_cols = ['temp']
X = bikes[feature_cols]
y = bikes.total

# 导入，实例化，拟合
from sklearn.linear_model import LinearRegression
linreg = LinearRegression()
linreg.fit(X, y)

#打印系数
print linreg.intercept_
print linreg.coef_

# 解释**拦截**（$ \ beta_0 $）：
#  - 当$ x $ = 0时，它是$ y $的值。
#  - 因此，当温度为0摄氏度时，这是预计的租金数量。
#  - **注意：**解释截取并不总是有意义的。 （为什么？）
# 解读**“临时”系数**（$ \ beta_1 $）：
#  - $ y $的变化除以$ x $的变化，或“斜率”。
#  - 因此，1摄氏度的温度升高**与9.17辆自行车的租金增加相关**。
#  - 这不是一个因果关系的陈述。
#  - $ \ beta_1 $如果温度上升与租金下降**相关，则**为负**。

# ＃##使用模型进行预测
# ＃如果温度是25摄氏度，我们可以预测多少辆自行车出租？

# 手动计算预测
linreg.intercept_ + linreg.coef_*25

# 使用预测方法
linreg.predict(25)

# ##功能的规模是否重要？
# 假设温度是用华氏度来衡量的，而不是摄氏度。 这将如何影响模型？
# 为华氏温度创建一个新的列
bikes['temp_F'] = bikes.temp * 1.8 + 32
bikes.head()

# Seaborn散点图与回归线
sns.lmplot(x='temp_F', y='total', data=bikes, aspect=1.5, scatter_kws={'alpha':0.2})

# 创建X和Y.
feature_cols = ['temp_F']
X = bikes[feature_cols]
y = bikes.total

# 实例化和拟合
linreg = LinearRegression()
linreg.fit(X, y)

# 打印系数
print linreg.intercept_
print linreg.coef_

# 将华氏温度转换成华氏25度
25 * 1.8 + 32

# 预测租金为华氏77度
linreg.predict(77)

# **结论：**线性回归模型的特征尺度是**不相关**。 在改变比例尺的时候，我们只是简单的改变**系数的**解释。
# 删除temp_F列
bikes.drop('temp_F', axis=1, inplace=True)

# 探索更多功能
feature_cols = ['temp', 'season', 'weather', 'humidity']

#在Seaborn多个散点图
sns.pairplot(bikes, x_vars=feature_cols, y_vars='total', kind='reg')

# 熊猫中的多个散点图
fig, axs = plt.subplots(1, len(feature_cols), sharey=True)
for index, feature in enumerate(feature_cols):
    bikes.plot(kind='scatter', x=feature, y='total', ax=axs[index], figsize=(16, 3))

# 你有没有看到你没有想到的事情？

# 季节和月份的交叉列表
pd.crosstab(bikes.season, bikes.index.month)

# 方块租金，按季节分组
bikes.boxplot(column='total', by='season')

# 值得注意的是：
#  - 行不能捕获非线性关系。
#  - 冬天有更多的租金比春天（？）
bikes.total.plot()

# 这告诉我们什么？
# 冬天比春天有更多的租金，但仅仅是因为系统正在经历**整体增长**而冬季月份恰好在春季之后。
# 相关矩阵（范围从1到-1）
bikes.corr()#计算列的成对相关性，不包括NA /空值

# 使用热图可视化Seaborn中的相关矩阵
sns.heatmap(bikes.corr())

# 你注意到了什么关系？
# ##为模型添加更多功能
# 创建一个功能列表

feature_cols = ['temp', 'season', 'weather', 'humidity']

# 创建X和Y.
X = bikes[feature_cols]
y = bikes.total

# 实例化和拟合
linreg = LinearRegression()
linreg.fit(X, y)

# 打印系数
print linreg.intercept_
print linreg.coef_

# 将要素名称与系数配对
zip(feature_cols, linreg.coef_)

# ＃解读系数：
#  - 固定所有其他功能，**温度增加1个单位**与**租赁增加7.86个自行车相关**。
#  - 固定所有其他功能，在**季节增加1个单位**与**租赁增加22.5个自行车相关**。
#  - 固定所有其他功能，**天气增加1个单位**与**租金增加6.67个自行车相关**。
#  - 固定所有其他功能，**湿度增加1个单位**与**减少3.12个自行车** **相关。

# 有什么看起来不正确？
###功能选择
# 我们如何选择模型中包含哪些功能？我们将要使用** train / test split **（最终**交叉验证**）。

# 为什么不使用** p值**或** R平方**来进行特征选择？
#  - 线性模型依赖**很多假设**（如特征是独立的），如果违反这些假设，则p值和R平方不太可靠。训练/测试分割依赖于较少的假设。
#  - 与响应无关的特征仍然可以具有**显着的p值**。
#  - 向你的模型添加与响应无关的特征将总是**增加R平方值**，而调整后的R平方不能充分解释这一点。

# 对于我们的泛化目标，P值和R平方是**代理**，而列车/测试分割和交叉验证尝试**直接估计**模型将如何推广到样本外数据。
#  更普遍：
#  - 有不同的方法可以用来解决任何给定的数据科学问题，本课程遵循**机器学习方法**。
#  - 本课程重点讨论**通用方法**，可以应用于任何模型，而不是模型特定的方法。

###回归问题的评估指标
# 分类问题的评估指标，如**准确性**，对回归问题没有用处。我们需要评估指标来比较**连续值**。
# 下面是回归问题的三个常见评估指标：
# **平均绝对误差**（MAE）是误差绝对值的平均值：
# $$ \ frac 1n \ sum_ {i = 1} ^ n | y_i- \ hat {y} _i | $$
# **均方误差**（MSE）是平方误差的均值：
# $$ \ frac 1n \ sum_ {i = 1} ^ n（y_i- \ hat {y} _i）^ 2 $$
# **均方根误差**（RMSE）是平方误差平均值的平方根：

# ＃$$ \ sqrt {\ frac 1n \ sum_ {i = 1} ^ n（y_i- \ hat {y} _i）^ 2} $$
#
# ＃示例真实和预测的响应值
true = [10, 7, 5, 5]
pred = [8, 6, 5, 10]

# 手工计算这些指标！
from sklearn import metrics
import numpy as np
print 'MAE:', metrics.mean_absolute_error(true, pred)
print 'MSE:', metrics.mean_squared_error(true, pred)
print 'RMSE:', np.sqrt(metrics.mean_squared_error(true, pred))

# 比较这些指标：
#  - ** MAE **是最容易理解的，因为这是平均误差。
#  - ** MSE **比MAE更受欢迎，因为MSE“惩罚”更大的错误，这在现实世界中往往是有用的。
#  - ** RMSE **比MSE更受欢迎，因为RMSE可以用“y”单位来解释。
# 所有这些都是**损失函数**，因为我们想要最小化它们。
# 这是另一个示例，演示MSE / RMSE如何惩罚更大的错误：
# 与上面相同的真值

true = [10, 7, 5, 5]

# 新的一组预测值
pred = [10, 7, 5, 13]

# MAE和以前一样
print 'MAE:', metrics.mean_absolute_error(true, pred)

# MSE和RMSE比以前更大
print 'MSE:', metrics.mean_squared_error(true, pred)
print 'RMSE:', np.sqrt(metrics.mean_squared_error(true, pred))

# ##将模型与训练/测试分离和RMSE进行比较

from sklearn.cross_validation import train_test_split

# 定义一个接受功能列表并返回测试RMSE的函数
def train_test_rmse(feature_cols):
    X = bikes[feature_cols]
    y = bikes.total
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=123)
    linreg = LinearRegression()
    linreg.fit(X_train, y_train)
    y_pred = linreg.predict(X_test)
    return np.sqrt(metrics.mean_squared_error(y_test, y_pred))

# 比较不同的功能集
print train_test_rmse(['temp', 'season', 'weather', 'humidity'])
print train_test_rmse(['temp', 'season', 'weather'])
print train_test_rmse(['temp', 'season', 'humidity'])

# 使用这些功能是不允许的！
print train_test_rmse(['casual', 'registered'])

# ##比较测试RMSE和空RMSE
# 无效RMSE是可以通过总是预测平均响应值**来实现的RMSE。 这是您可能想要衡量您的回归模型的基准。
# 把X和Y分成训练和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=123)

# 创建一个与y_test形状相同的NumPy数组
y_null = np.zeros_like(y_test, dtype=float)

# 用y_test的平均值填充数组
y_null.fill(y_test.mean())
y_null

# 计算空值RMSE
np.sqrt(metrics.mean_squared_error(y_test, y_null))

###处理分类特征
# scikit-learn期望所有功能都是数字。 那么我们如何在模型中包含分类特征呢？
#  - **有序类别：**将它们转换为合理的数值（例如：small = 1，medium = 2，large = 3）
#  - **无序的类别：**使用虚拟编码（0/1）
#数据集中的分类特征是什么？
#  - **有序类别：**天气（已经用合理的数值编码）
#  - **无序类别：**季节（需要虚拟编码），节假日（已经虚拟编码），工作日（已经虚拟编码）
# 对于季节，我们不能简单地将编码保留为1 =春天，2 =夏天，3 =秋天，4 =冬天，因为那意味着**有序的关系**。
# 相反，我们创建**多个虚拟变量：**

# 创建虚拟变量
season_dummies = pd.get_dummies(bikes.season, prefix='season')

# 打印5个随机行
season_dummies.sample(n=5, random_state=1)

# 但是，我们实际上只需要**三个虚拟变量（而不是四个）**，因此我们将放弃第一个虚拟变量。
# 为什么？ 因为三个假人捕捉了关于季节特征的所有“信息”，并且隐含地将春季（季节1）定义为**基线水平：**
# ＃放下第一列
season_dummies.drop(season_dummies.columns[0], axis=1, inplace=True)

#打印5个随机行
season_dummies.sample(n=5, random_state=1)

# 一般来说，如果您有** k个可能值**的分类特征，则可以创建** k-1个虚拟变量**。
#如果这是令人困惑的，那么考虑一下为什么我们只需要一个虚拟变量来放假，而不是两个虚拟变量（holiday_yes和holiday_no）。
# ＃连接原始DataFrame和虚拟DataFrame（axis = 0表示行，axis = 1表示列）
bikes = pd.concat([bikes, season_dummies], axis=1)

# 打印5个随机行
bikes.sample(n=5, random_state=1)

# 在模型中包含季节的虚拟变量
feature_cols = ['temp', 'season_2', 'season_3', 'season_4', 'humidity']
X = bikes[feature_cols]
y = bikes.total
linreg = LinearRegression()
linreg.fit(X, y)
zip(feature_cols, linreg.coef_)

# 我们如何解读季节系数？ 他们是**相对于基线（春季）**测量的：
#  - 保持所有其他功能固定，**夏季**与**租金减少3.39辆自行车**相比，春天。
#  - 保持所有其他功能固定，**下降**与租金减少41.7辆自行车**相比，春天。
#  - 保持所有其他功能固定，**冬季**与64.4辆自行车**相比，春季出租增加**。
# 如果我们改变哪个赛季被定义为基线，这会有影响吗？
# - 不，它只会改变我们**系数的解释。
# **重要：**虚拟编码与所有机器学习模型相关，而不仅仅是线性回归模型。
# 比较原始季节变量与虚拟变量
print train_test_rmse(['temp', 'season', 'humidity'])
print train_test_rmse(['temp', 'season_2', 'season_3', 'season_4', 'humidity'])

# ##特征工程
#看看你是否可以创建以下功能：
#  - **小时：**作为单个数字功能（0到23）
#  - **小时：**作为分类特征（使用23个虚拟变量）
# - **白天：**作为单个分类特征（白天=上午7点到晚上8点，否则白天= 0）
# 然后，尝试使用`train_test_rmse`三个功能中的每一个（用它自己）来看哪一个功能最好！
# 小时作为数字功能
bikes['hour'] = bikes.index.hour

# 小时作为分类特征
hour_dummies = pd.get_dummies(bikes.hour, prefix='hour')
hour_dummies.drop(hour_dummies.columns[0], axis=1, inplace=True)
bikes = pd.concat([bikes, hour_dummies], axis=1)

# 白天作为一个明确的特征
bikes['daytime'] = ((bikes.hour > 6) & (bikes.hour < 21)).astype(int)

print train_test_rmse(['hour'])
print train_test_rmse(bikes.columns[bikes.columns.str.startswith('hour_')])
print train_test_rmse(['daytime'])

# ##比较线性回归与其他模型
# 线性回归的优点：
#  - 简单解释
#  - 高度可解释的
#  - 模型训练和预测是快速的
#  - 不需要调整（不包括正规化）
#  - 功能不需要缩放
#  - 可以在少量的观察中表现良好
#  - 完全了解

# 线性回归的缺点：
#  - 设定功能和响应之间的线性关系
#  - 由于高偏见，表现（通常）与最好的监督式学习方法没有竞争力

第二课 用Yelp投票线性回归作业

# 介绍
# 此作业使用Kaggle的[Yelp Business Rating Prediction]（https://www.kaggle.com/c/yelp-recsys-2013）竞赛中的一小部分数据。

#  - 此数据集中的每个观察结果都是特定用户对特定业务的评论。
#  - “星号”列是审阅者分配给企业的星号（1到5）。 （更高的星星更好。）换句话说，这是撰写评论的人的业务评级。
#  - “酷”栏是本评论从其他Yelp用户收到的“酷”票数。
# 所有的评论都是从0票“冷静”开始的，评论可以获得多少“酷”的票数是没有限制的。
# 换句话说，这是对评论本身的评价，而不是对评论的评价。
#  - “有用”和“有趣”列与“酷”栏相似。

# 将`yelp.csv`读入DataFrame。
# 使用相对路径访问yelp.csv

import pandas as pd
yelp = pd.read_csv('yelp.csv')
yelp.head(1)

###任务1（奖金）
# 忽略`yelp.csv`文件，并从`yelp.json`自己构建这个DataFrame。
# 这包括将数据读入Python，解码JSON，将其转换为DataFrame，并为每个投票类型添加单独的列。
# 从yelp.json读取数据到行列表中
# 使用json.loads（）将每行解码成字典
import json
with open('yelp.json', 'rU') as f:
    data = [json.loads(row) for row in f]

# 显示第一个评论
data[0]

# 将字典列表转换为DataFrame
yelp = pd.DataFrame(data)
yelp.head(1)

#添加DataFrame列以获得很酷，有用和有趣的内容
yelp['cool'] = [row['votes']['cool'] for row in data]
yelp['useful'] = [row['votes']['useful'] for row in data]
yelp['funny'] = [row['votes']['funny'] for row in data]

#放弃选票栏
yelp.drop('votes', axis=1, inplace=True)
yelp.head(1)

# ＃##任务2
#探索每种投票类型（酷/有用/有趣）和星星数量之间的关系。
#把星星当成一个分类变量，寻找组间的差异
yelp.groupby('stars').mean()

#相关矩阵
import seaborn as sns
sns.heatmap(yelp.corr())

# 多个散点图
sns.pairplot(yelp, x_vars=['cool', 'useful', 'funny'], y_vars='stars', size=6, aspect=0.7, kind='reg')

# ## Task 3
#
# 定义酷/有用/有趣的功能，明星作为响应。

feature_cols = ['cool', 'useful', 'funny']
X = yelp[feature_cols]
y = yelp.stars

# ## Task 4
#
# 拟合线性回归模型并解释系数。 这些系数对你来说有直觉意义吗？ 浏览Yelp网站，查看是否检测到类似的趋势。
from sklearn.linear_model import LinearRegression
linreg = LinearRegression()
linreg.fit(X, y)
zip(feature_cols, linreg.coef_)

# ## Task 5
#
# 通过将模型分解为训练集和测试集并计算RMSE来评估模型。 RMSE是否对你有直觉意义？
from sklearn.cross_validation import train_test_split
from sklearn import metrics
import numpy as np

# 定义一个接受功能列表并返回测试RMSE的函数

# 做回归分析，常用的误差主要有均方误差根（RMSE）和R-平方（R2）。
# RMSE是预测值与真实值的误差平方根的均值
# R2方法是将预测值跟只使用均值的情况下相比，看能好多少

def train_test_rmse(feature_cols):
    X = yelp[feature_cols]
    y = yelp.stars
    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)
    linreg = LinearRegression()
    linreg.fit(X_train, y_train)
    y_pred = linreg.predict(X_test)
    return np.sqrt(metrics.mean_squared_error(y_test, y_pred))

# 用所有三个特征计算RMSE
train_test_rmse(['cool', 'useful', 'funny'])

# ## Task 6
#
# 尝试删除一些功能，看看RMSE是否改善。

print train_test_rmse(['cool', 'useful'])
print train_test_rmse(['cool', 'funny'])
print train_test_rmse(['useful', 'funny'])

# ## Task 7 (Bonus)
#
# 思考一下可以从现有数据创建的一些新功能，这些功能可以预测响应。
# 找出如何在Pandas中创建这些功能，将它们添加到您的模型中，并查看RMSE是否得到改善。
# 新功能：评论长度（字符数）
yelp['length'] = yelp.text.apply(len)

# 新功能：评论是否包含“爱”或“讨厌”
yelp['love'] = yelp.text.str.contains('love', case=False).astype(int)
yelp['hate'] = yelp.text.str.contains('hate', case=False).astype(int)

# 为模型添加新功能并计算RMSE
train_test_rmse(['cool', 'useful', 'funny', 'length', 'love', 'hate'])

# ## Task 8 (Bonus)
#
# 把你最好的RMSE和测试集合中的RMSE比较为“空模型”，这是忽略所有特征的模型，并简单地预测测试集合中的平均响应值。

#分割数据（在函数之外）
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

# 创建一个与y_test形状相同的NumPy数组
y_null = np.zeros_like(y_test, dtype=float)
# 用y_test的均值填充数组
y_null.fill(y_test.mean())
# 计算零均方根误差
print np.sqrt(metrics.mean_squared_error(y_test, y_null))

# ## Task 9 (Bonus)
#
#不要把它看作是一个回归问题，而应该把它看作一个分类问题，看看用KNN可以达到什么样的测试精度。
#导入并实例化KNN
from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=50)

# 分类模型将自动将响应值（1/2/3/4/5）视为无序类别knn.fit(X_train, y_train)
knn.fit(X_train, y_train)
y_pred_class = knn.predict(X_test)
print metrics.accuracy_score(y_test, y_pred_class)

# ## Task 10 (Bonus)
# 找出如何使用线性回归进行分类，并将其分类精度与KNN的准确度进行比较。

# 使用线性回归进行连续预测
linreg = LinearRegression()
linreg.fit(X_train, y_train)
y_pred = linreg.predict(X_test)

# 将其预测整理为最接近的整数
y_pred_class = y_pred.round()

# 计算舍入预测的分类准确度
print metrics.accuracy_score(y_test, y_pred_class)

第三课：指数函数和对数

import math
import numpy as np

# ## 指数函数

# 什么是** e **？ 这只是一个数字（被称为欧拉数）：

math.e

#** e **是一个重要的数字，因为它是所有持续增长的流程共享的基础增长率。
#例如，如果我有** 10美元**，而且在1年内增长100％（持续复利），我会以** 10 \ * e ^ 1美元结束**：

# 1年增长100％
10 * np.exp(1)
# 100％增长2年
10 * np.exp(2)

#注意：当e上升到一个幂时，它被称为**指数函数**。
# 从技术上讲，任何数字都可以作为基数，它仍然被称为**指数函数**（如2 ^ 5）。
# 但是在我们的情况下，指数函数的基础被假定为e。
#无论如何，如果我只有20％的增长，而不是100％的增长呢？

# 1年增长20％
10 * np.exp(0.20)
# 2年增长20％
10 * np.exp(0.20 * 2)

###对数
#（自然对数）**是什么？ 它给你所需的时间达到一定的增长水平。
# 例如，如果我想增长2.718倍，则需要1个单位的时间（假设增长率为100％）：
#需要增加1个单位到2.718个单位

np.log(2.718)

#如果我要增长7.389倍，那么我需要2个单位的时间：
# 需要增加1个单位到7.389个单位
np.log(7.389)

# 如果我想要增长1倍，它会花费我0个单位的时间：
# 需要将1个单位增加到1个单位的时间
np.log(1)

#如果我想要增长0.5倍，它将花费我-0.693个单位的时间（这就像是在回顾时间）：
#需要增加1个单位到0.5个单位
np.log(0.5)

##连接概念
#正如你所看到的那样，指数函数和自然对数是**反转**：

np.log(np.exp(5))
np.exp(np.log(5))

第四课、 Logistic回归

#--logistic回归又称logistic回归分析，是一种广义的线性回归分析模型，常用于数据挖掘，疾病自动诊断，经济预测等领域

# ## Part 1: 预测连续响应

# 玻璃鉴定数据集
import pandas as pd
url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/glass/glass.data';;
col_names = ['id','ri','na','mg','al','si','k','ca','ba','fe','glass_type']
glass = pd.read_csv(url, names=col_names, index_col='id')
glass.sort_values('al', inplace=True)
glass.head()

# **问题：**假装我们要预测** ri **，而我们唯一的特征是** al **。 我们怎么能使用机器学习呢？
# **答案：**我们可以将其作为一个回归问题，并使用** al **作为唯一特征和** ri **作为响应的线性回归模型。

# **问题：**我们如何将**这个模型可视化？
# **答案：**用x轴上的** al **和y轴上的** ri **创建一个散点图，并画出最适合的线。

import seaborn as sns
import matplotlib.pyplot as plt
sns.set(font_scale=1.5)

sns.lmplot(x='al', y='ri', data=glass, ci=None)

#**问题：**我们如何在不使用Seaborn的情况下画出这个情节？

# 散点图使用熊猫
glass.plot(kind='scatter', x='al', y='ri')

# 使用Matplotlib的等效散点图
plt.scatter(glass.al, glass.ri)
plt.xlabel('al')
plt.ylabel('ri')

# 拟合线性回归模型
from sklearn.linear_model import LinearRegression
linreg = LinearRegression()
feature_cols = ['al']
X = glass[feature_cols]
y = glass.ri
linreg.fit(X, y)

# 预测X的所有值
glass['ri_pred'] = linreg.predict(X)
glass['ri_pred'].head()

# 绘制一条线连接的预测
plt.plot(glass.al, glass.ri_pred, color='red')
plt.xlabel('al')
plt.ylabel('Predicted ri')

# 把这些情节放在一起
plt.scatter(glass.al, glass.ri)
plt.plot(glass.al, glass.ri_pred, color='red')
plt.xlabel('al')
plt.ylabel('ri')

###进修：解释线性回归系数

#线性回归方程：$ y = \ beta_0 + \ beta_1x $

#使用公式计算al = 2的预测

linreg.intercept_ + linreg.coef_ * 2

# 使用预测方法计算al = 2的预测
linreg.predict(2)

# 检查系数
zip(feature_cols, linreg.coef_)

#**解释：**'al'增加1个单位与'ri'减少0.0025个单位相关。

#将al加1（使得al = 3）将ri减小0.0025
1.51699012 - 0.0024776063874696243

# 使用预测方法计算对于al = 3的预测
linreg.predict(3)

##第2部分：预测分类响应

#检查glass_type
glass.glass_type.value_counts().sort_index()

#1，2，3是窗户玻璃
#5,6,7型是家用玻璃
glass['household'] = glass.glass_type.map({1:0, 2:0, 3:0, 5:1, 6:1, 7:1})
glass.head()

# ＃让我们改变我们的任务，以便我们使用** al **预测**家庭**。 让我们通过可视化的关系找出如何做到这一点：
plt.scatter(glass.al, glass.household)
plt.xlabel('al')
plt.ylabel('household')

# 让我们画一条**回归线**，就像我们之前做的那样：

# 拟合线性回归模型并存储预测
feature_cols = ['al']
X = glass[feature_cols]
y = glass.household
linreg.fit(X, y)
glass['household_pred'] = linreg.predict(X)

#包含回归线的散点图
plt.scatter(glass.al, glass.household)
plt.plot(glass.al, glass.household_pred, color='red')
plt.xlabel('al')
plt.ylabel('household')

#如果** al = 3 **，我们预测哪一类家庭？****1
#如果** al = 1.5 **，我们预测什么样的家庭类？**0**
#我们预测了al的** lower **值的0类，以及al的** higher **值的1个类。
#我们的截止值是多少？ 在** al = 2 **附近，因为这是线性回归线跨越预测0级和1级之间的中点的地方。

#因此，如果** household_pred> = 0.5 **，我们预测** 1 **类，否则我们预测** 0 **类。

#理解np.where
import numpy as np
nums = np.array([5, 15, 8])

# 如果条件为True，则np.where返回第一个值，如果条件为False，则返回第二个值
np.where(nums > 10, 'big', 'small')

# 将household_pred转换为1或0
glass['household_pred_class'] = np.where(glass.household_pred >= 0.5, 1, 0)
glass.head()

# 绘制类别预测
plt.scatter(glass.al, glass.household)
plt.plot(glass.al, glass.household_pred_class, color='red')
plt.xlabel('al')
plt.ylabel('household')

###第三部分：使用Logistic回归
#逻辑回归可以做我们刚刚做的事情：

#适合逻辑回归模型并存储类别预测
from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression(C=1e9)
feature_cols = ['al']
X = glass[feature_cols]
y = glass.household
logreg.fit(X, y)
glass['household_pred_class'] = logreg.predict(X)

# 绘制类别预测
plt.scatter(glass.al, glass.household)
plt.plot(glass.al, glass.household_pred_class, color='red')
plt.xlabel('al')
plt.ylabel('household')

#如果我们想要**预测概率**，而不是仅仅**类预测**，那么理解我们在给定预测中有多自信？

#存储第1类的预测概率
glass['household_pred_prob'] = logreg.predict_proba(X)[:, 1]

# 绘制预测概率
plt.scatter(glass.al, glass.household)
plt.plot(glass.al, glass.household_pred_prob, color='red')
plt.xlabel('al')
plt.ylabel('household')

# 检查一些示例预测
print logreg.predict_proba(1)
print logreg.predict_proba(2)
print logreg.predict_proba(3)

###第4部分：概率，赔率，电子，日志，对数

# 创建一个概率与赔率的表
table = pd.DataFrame({'probability':[0.1, 0.2, 0.25, 0.5, 0.6, 0.8, 0.9]})
table['odds'] = table.probability/(1 - table.probability)
table

#什么是** e **？ 这是所有持续增长的进程共享的基本增长率：

#指数函数：e ^ 1
np.exp(1)

#什么是**（自然）日志**？ 它给你所需的时间达到一定的增长水平：

#需要增加1个单位到2.718个单位
np.log(2.718)

# 这也是指数函数的**反**。

np.log(np.exp(5))

# 添加日志赔率到表中
table['logodds'] = np.log(table.odds)
table

# ## ＃##第五部分：什么是Logistic回归？

# **线性回归：**连续响应被建模为这些特征的线性组合：
#
# $$y = \beta_0 + \beta_1x$$
#
# 逻辑回归：**分类回答为“真”（1）的对数可能性被建模为这些特征的线性组合#
# $$\log \left({p\over 1-p}\right) = \beta_0 + \beta_1x$$
#
# 这被称为** logit函数**。

###第六部分：解释Logistic回归系数
#再次绘制预测概率
plt.scatter(glass.al, glass.household)
plt.plot(glass.al, glass.household_pred_prob, color='red')
plt.xlabel('al')
plt.ylabel('household')

# 使用equ计算al = 2的预测对数
logodds = logreg.intercept_ + logreg.coef_[0] * 2
logodds

# 将对数赔率转换为赔率
odds = np.exp(logodds)
odds

# 赔率转换为概率
prob = odds/(1 + odds)
prob

# 使用predict_proba方法计算al = 2的预测概率
logreg.predict_proba(2)[:, 1]

#检查al的系数
zip(feature_cols, logreg.coef_[0])

#**解释：**“1”增加1个单位，与“家庭”的对数增加4.18个单位相关。
#增加1（这样al = 3）使对数增加4.18
logodds = 0.64722323 + 4.1804038614510901
odds = np.exp(logodds)
prob = odds/(1 + odds)
prob

# 使用predict_proba方法计算al = 3的预测概率
logreg.predict_proba(3)[:, 1]

#**底线：**正系数增加响应的对数（从而增加概率），负系数减少响应的对数（从而降低概率）。
#检查拦截
logreg.intercept_

#**解释：**对于“al”值为0，“家庭”的对数为-7.71。
#把log-odds转换成概率

logodds = logreg.intercept_
odds = np.exp(logodds)
prob = odds/(1 + odds)
prob

###第7部分：使用具有分类特征的Logistic回归

#逻辑回归仍然可以与**分类特征**一起使用。 让我们看看是什么样的：

#创建一个分类功能
glass['high_ba'] = np.where(glass.ba > 0.5, 1, 0)

#使用Seaborn绘制逻辑曲线：

#原创（连续）功能
sns.lmplot(x='ba', y='household', data=glass, ci=None, logistic=True)

# 分类特征
sns.lmplot(x='high_ba', y='household', data=glass, ci=None, logistic=True)

# 分类功能，添加了抖动
sns.lmplot(x='high_ba', y='household', data=glass, ci=None, logistic=True, x_jitter=0.05, y_jitter=0.05)

# 拟合逻辑回归模型
feature_cols = ['high_ba']
X = glass[feature_cols]
y = glass.household
logreg.fit(X, y)

# 检查high_ba的系数
zip(feature_cols, logreg.coef_[0])

###第八部分：比较Logistic回归与其他模型
# 逻辑回归的优点：
#  - 高度可解释的（如果你还记得的话）
#  - 模型训练和预测是快速的
#  - 不需要调整（不包括正规化）
#  - 功能不需要缩放
#  - 可以在少量的观察中表现良好
#  - 输出经过良好校准的预测概率

# 逻辑回归的缺点：
#  - 设定特征和响应的对数之间的线性关系
#  - 表现（通常）与最好的监督式学习方法没有竞争力
#  - 无法自动学习功能交互

第五课与泰坦尼克号数据进行逻辑回归练习

# # 用泰坦尼克数据进行逻辑回归练习

# ## Step 1: 将数据读入Pasnda

import pandas as pd
titanic = pd.read_csv('titanic.csv', index_col='PassengerId')
titanic.head()

# ## Step 2: 创建X和Y
#
# 定义** Pclass **和** Parch **作为特征，** Survived **作为响应。

feature_cols = ['Pclass', 'Parch']
X = titanic[feature_cols]
y = titanic.Survived

# ## Step 3: 将数据分解成训练和测试集。

from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)

# ## Step 4: 拟合逻辑回归模型并检查系数
#
# 确认系数具有直觉意义。

from sklearn.linear_model import LogisticRegression
logreg = LogisticRegression(C=1e9)
logreg.fit(X_train, y_train)
zip(feature_cols, logreg.coef_[0])

# ## Step 5: 对测试集进行预测并计算精度

# 类预测（不是预测概率）
y_pred_class = logreg.predict(X_test)

# 计算分类精度
from sklearn import metrics
print metrics.accuracy_score(y_test, y_pred_class)

# ## Step 6: 将您的测试准确性与零准确度进行比较

# 这个工程不管类的数量
y_test.value_counts().head(1) / len(y_test)

# 这只适用于编码为0/1的二进制分类问题
max(y_test.mean(), 1 - y_test.mean())

#泰坦尼克号预测的混淆矩阵＃

#打印混淆矩阵
print metrics.confusion_matrix(y_test, y_pred_class)

#保存混淆矩阵并切片成四块
confusion = metrics.confusion_matrix(y_test, y_pred_class)
TP = confusion[1][1]
TN = confusion[0][0]
FP = confusion[0][1]
FN = confusion[1][0]

print 'True Positives:', TP
print 'True Negatives:', TN
print 'False Positives:', FP
print 'False Negatives:', FN

# 计算灵敏度
print TP / float(TP + FN)
print 44 / float(44 + 51)

#计算特异性
print TN / float(TN + FP)
print 105 / float(105 + 23)

# 存储预测的概率
y_pred_prob = logreg.predict_proba(X_test)[:, 1]

# 预测概率的直方图
import matplotlib.pyplot as plt
plt.hist(y_pred_prob)
plt.xlim(0, 1)
plt.xlabel('Predicted probability of survival')
plt.ylabel('Frequency')

# 通过降低预测生存的阈值来提高灵敏度
import numpy as np
y_pred_class = np.where(y_pred_prob > 0.3, 1, 0)

# 旧混淆矩阵
print confusion

# 新的混淆矩阵
print metrics.confusion_matrix(y_test, y_pred_class)

#新的灵敏度（比以前更高）
print 63 / float(63 + 32)

# 新的特异性（低于以前）
print 72 / float(72 + 56)

第一课：从电影网站的API获取数据

# 将IMDb数据读入DataFrame中：我们需要一年的专栏！
import pandas as pd
import requests
from time import sleep
from bs4 import BeautifulSoup

movies = pd.read_csv('imdb_1000.csv')
movies.head()

# 使用请求库来与URL进行交互
r = requests.get('http://www.omdbapi.com/?apikey=19a14323&t=the shawshank redemption&r=json&type=movie')

# 检查状态：200表示成功，4xx表示错误
r.status_code

# 查看原始响应文本
r.text

# 查看原始响应文本
r.json()

# 从字典中提取年
r.json()['Year']

# 如果电影名称不被识别会怎么样？?
r = requests.get('http://www.omdbapi.com/?apikey=19a14323&t=blahblahblah&r=json&type=movie';;)
r.status_code
r.json()

# 定义一个返回年份的函数
def get_movie_year(title):
    r = requests.get('http://www.omdbapi.com/?apikey=19a14323&t=';; + title + '&r=json&type=movie')
    info = r.json()
    if info['Response'] == 'True':
        return int(info['Year'])
    else:
        return None

# 测试功能
get_movie_year('The Shawshank Redemption')
get_movie_year('blahblahblah')

# 创建一个较小的DataFrame进行测试
top_movies = movies.head().copy()

# 写一个for循环来建立一个年份列表
years = []
for title in top_movies.title:
    years.append(get_movie_year(title))
    sleep(1)

# 检查DataFrame和年份列表的长度是否相同
assert(len(top_movies) == len(years))

# 将该列表保存为新列
top_movies['year'] = years

'''
Bonus 更新DataFrame作为循环的一部分
'''

# 枚举允许您在迭代时访问项目位置
letters = ['a', 'b', 'c']
for index, letter in enumerate(letters):
    print index, letter

#DataFrames的iterrows方法是类似的
for index, row in top_movies.iterrows():
    print index, row.title

# 创建一个新的列并设置一个默认值
movies['year'] = -1

# loc方法允许您通过'label'访问DataFrame元素
movies.loc[0, 'year'] = 1994

#写一个for循环来更新前三部电影的年份
for index, row in movies.iterrows():
    if index < 3:
        movies.loc[index, 'year'] = get_movie_year(row.title)
        sleep(1)
    else:
        break

第二课、从网站上爬取网页，并解析存入pandas

# 阅读网页的HTML代码并保存为字符串
with open('example.html', 'rU') as f:
    html = f.read()

# 将HTML转换为结构化的Soup对象

b = BeautifulSoup(html)

# 打印出对象
print b
print b.prettify()

    # 'find'方法返回第一个匹配的标签（以及其中的所有内容）
b.find(name='body')
b.find(name='h1')

# 标签允许您访问“内部文本”
b.find(name='h1').text

# 标签还允许您访问其属性
b.find(name='h1')['id']

# “find_all”方法对于查找所有匹配的标签非常有用
b.find(name='p')        # 返回一个标签
b.find_all(name='p')    # 返回一个ResultSet（就像一个标签列表）

# ResultSets可以像列表一样切片
len(b.find_all(name='p'))
b.find_all(name='p')[0]
b.find_all(name='p')[0].text
b.find_all(name='p')[0]['id']

#迭代ResultSet
results = b.find_all(name='p')
for tag in results:
    print tag.text

# 限制标签属性的搜索
b.find(name='p', attrs={'id':'scraping'})
b.find_all(name='p', attrs={'class':'topic'})
b.find_all(attrs={'class':'topic'})

#限制搜索到特定的部分
b.find_all(name='li')
b.find(name='ul', attrs={'id':'scraping'}).find_all(name='li')

'''
练习一
'''

# 找到'h2'标签然后打印它的文本
b.find(name='h2').text

# 找到“feedback”的“id”值的“p”标签，然后打印其文本
b.find(name='p', attrs={'id':'feedback'}).text

# 找到第一个“p”标签，然后打印“id”属性的值
b.find(name='p')['id']

# 打印所有四个资源的文本
results = b.find_all(name='li')
for tag in results:
    print tag.text

# 只打印API资源的文本
results = b.find(name='ul', attrs={'id':'api'}).find_all(name='li')
for tag in results:
    print tag.text

# 从Shawshank Redemption页面获取HTML

r = requests.get('http://www.imdb.com/title/tt0111161/')

# 将HTML转换成soup
b = BeautifulSoup(r.text)
print b

# 如果您有编码错误，请运行此代码
import sys
reload(sys)
sys.setdefaultencoding('utf8')

b.find_all(name='span', attrs={'class':'itemprop', 'itemprop':'name'})    #结果太多了
b.find(name='span', attrs={'class':'itemprop', 'itemprop':'name'}).text   # 获得导演

b.find(name='div',attrs={'class':'title_wrapper'}).find(name='h1').text# 得到标题

# 将星级（作为一个浮点数）

float(b.find(name='span', attrs={'itemprop':'ratingValue'}).text)

'''
练习二
'''

# 得到描述
b.find(name='div', attrs={'itemprop':'description'}).text.strip()

# 获得内容评级
b.find(name='meta', attrs={'itemprop':'contentRating'})['content']

#获得分钟的持续时间（作为一个整数）
str(b.find(name='time', attrs={'itemprop':'duration'})['datetime'][2:-1])

'''
定义一个接受IMDb ID并返回一个字典的函数
电影信息：标题，星级，说明，content_rating，持续时间
'''

# 定义一个接受IMDb ID并返回电影信息字典的函数
def get_movie_info(imdb_id):
    r = requests.get('http://www.imdb.com/title/';; + imdb_id + '/')
    b = BeautifulSoup(r.text)
    info = {}
    info['title'] = b.find(name='div',attrs={'class':'title_wrapper'}).find(name='h1').text
    info['star_rating'] = float(b.find(name='span', attrs={'itemprop':'ratingValue'}).text)
    info['description'] = b.find(name='div', attrs={'itemprop':'description'}).text.strip()
    info['content_rating'] = b.find(name='meta', attrs={'itemprop':'contentRating'})['content']
    info['duration'] = str(b.find(name='time', attrs={'itemprop':'duration'})['datetime'][2:-1])
    return info

# 测试方法
get_movie_info('tt0109830')

# 打开ID文件（每行一个ID），并将ID存储在一个列表中
imdb_ids = []
with open('imdb_ids.txt', 'rU') as f:
    imdb_ids = [row.strip() for row in f]

imdb_ids

# 获取每部电影的信息，并将结果存储在列表中
from time import sleep
movies = []
for imdb_id in imdb_ids:
    movies.append(get_movie_info(imdb_id))
    sleep(1)

movies

# 检查ID列表和电影列表长度是否相同
assert(len(imdb_ids) == len(movies))

# 将电影列表转换为DataFrame
import pandas as pd

pd.DataFrame(movies, index=imdb_ids)

第三课、探索偏差 - 方差权衡

import pandas as pd
import numpy as np
import seaborn as sns

# ## 大脑和体重

# 这是一个[数据集]（http://people.sc.fsu.edu/~jburkardt/datasets/regression/x01.txt）62种哺乳动物的身体和大脑的平均重量。
# 让我们把它读入熊猫，并快速浏览一下：

url = 'http://people.sc.fsu.edu/~jburkardt/datasets/regression/x01.txt';;
col_names = ['id', 'brain', 'body']
mammals = pd.read_table(url, sep='\s+', skiprows=33, names=col_names, index_col='id') #skiprowss 跳过33行
mammals.head()

mammals.describe() #统计信息

# 我们将重点关注体重小于200的一个较小的子集：

# 只保留体重小于200的排
mammals = mammals[mammals.body < 200]
mammals.shape

# 我们现在假装现在只有51种哺乳动物。 换句话说，我们假装这是每个已知哺乳动物物种**的整个脑和体重数据集。
# 让我们创建一个散点图
sns.lmplot(x='body', y='brain', data=mammals, ci=None, fit_reg=False)

# 哺乳动物的脑和体重之间似乎有一定的关系。

# ## 做一个预测

# 现在让我们假装发现了一种新的哺乳动物物种**。
#  我们测量我们可以找到的这个物种的每个成员的体重，并计算**的平均体重。
# 我们要**预测这个物种的平均脑重**（而不是直接测量）。 我们如何做到这一点？
sns.lmplot(x='body', y='brain', data=mammals, ci=None)

#我们画了一条直线，似乎最能捕捉到大脑和体重之间的关系。因此，我们可以预测，我们的新物种的脑重约45，因为这是x = 100时的近似y值。

#这被称为“线性模型”或“线性回归模型”，我们将在未来的课堂上学习。
##从样本中进行预测
# 早些时候，我说这个数据集包含了每一个已知的哺乳动物物种。这很方便，但是在现实世界中，你所拥有的仅仅是一个数据样本**。
# 一个更现实的情况是只有大脑和体重（比方说）51个已知哺乳动物的一半。
# 当发现新的哺乳动物（体重为100）时，我们仍然要对脑重做出准确的预测，但是这个任务可能会更困难，因为我们没有所有的数据理想喜欢有。
# 让我们模拟这种情况，将51个观测值分别分配给** 1或者2个宇宙**：

# 设置一个随机的种子重现性
np.random.seed(12345)

# 随机分配每个观察宇宙1或宇宙2
mammals['universe'] = np.random.randint(1, 3, len(mammals))
mammals.head()

# ＃**重要：**我们只住在两个宇宙中的一个。
# 这两个宇宙都有51种已知的哺乳动物物种，但是每个宇宙都知道不同物种的大脑和体重。
# ＃我们现在可以告诉Seaborn创建两个图，其中左边的图仅使用来自** universe 1 **的数据，右边的图只使用来自** universe 2 **的数据：

# col ='universe'按宇宙对数据进行分组，并创建两个独立的图
sns.lmplot(x='body', y='brain', data=mammals, ci=None, col='universe')

#这两条曲线看起来非常相似，尽管他们使用了单独的数据样本。 在这两种情况下，我们都可以预测脑重约45。

#通过将它们放在同一个图上，可以更容易地看到相似程度：

#hue ='universe'根据宇宙对数据进行子集并创建一个图
sns.lmplot(x='body', y='brain', data=mammals, ci=None, hue='universe')

# 这个练习的重点是什么？ 这是一个高偏见，低方差模型的视觉演示：
# 这是**高偏见**，因为它不适合数据特别好。
# 这是**低的方差**，因为它不会有太大的变化，取决于在那个宇宙中哪些观察结果是可用的。
# 让我们尝试一些完全不同的东西
# 低偏差，高方差**模型会是什么样子？ 让我们用一个八阶多项式来尝试多项式回归：

sns.lmplot(x='body', y='brain', data=mammals, ci=None, col='universe', order=8)

# ＃ - 这是**低偏见**因为模型匹配数据相当好！
# ＃ - 它是**高方差**，因为这些模型根据宇宙中哪些观察结果可用而有很大的不同。
# （对于100的体重，脑重预测将在一个宇宙中为40，在另一个宇宙中为0）
# ＃##我们能找到一个中间地带吗？
#
# ＃也许我们可以创建一个比**线性模型**偏差小的模型，**比八次多项式**更少**。
# ＃
# ＃我们来尝试一个二阶多项式：
sns.lmplot(x='body', y='brain', data=mammals, ci=None, col='universe', order=2)

第四课、K-最近的邻居和scikit-learn

# 将 iris data读入DataFrame
import pandas as pd
url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data';;
col_names = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'species']
iris = pd.read_csv(url, header=None, names=col_names)

iris.head()

import matplotlib.pyplot as plt

# 增加默认图形和字体大小以便于查看
plt.rcParams['figure.figsize'] = (6, 4)
plt.rcParams['font.size'] = 14

# 创建一个自定义颜色映射
from matplotlib.colors import ListedColormap
cmap_bold = ListedColormap(['#FF0000', '#00FF00', '#0000FF'])

# 将每个鸢尾花物种映射到一个数字
iris['species_num'] = iris.species.map({'Iris-setosa':0, 'Iris-versicolor':1, 'Iris-virginica':2})

# 根据物种创建花瓣长度与花瓣宽度和颜色的散点图
iris.plot(kind='scatter', x='petal_length', y='petal_width',
          c='species_num', colormap=cmap_bold)

# 创建一个萼片长度与萼片宽度和各种颜色的散点图
iris.plot(kind='scatter', x='sepal_length', y='sepal_width',
          c='species_num', colormap=cmap_bold)
iris.head()

# 将特征矩阵存储在“X”
feature_cols = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width']
X = iris[feature_cols]

# 替代方法来创建“X”
X = iris.drop(['species', 'species_num'], axis=1)
X = iris.loc[:, 'sepal_length':'petal_width']
X = iris.iloc[:, 0:4]

# store response vector in "y"
y = iris.species_num

# 检查X的类型
print type(X)
print type(X.values)

# 检查y的类型
print type(y)
print type(y.values)

# 检查X的形状（n =观察数量，p =特征数量）
print X.shape

# 检查y的形状（长度为n的单维）
print y.shape

# ## scikit-learn的4步建模模式

# **Step 1:** 导入你打算使用的类

from sklearn.neighbors import KNeighborsClassifier

# **步骤2：**“实例化”“估计器”
# - “Estimator”是scikit-learn的术语“模型”
# - “实例化”是指“创建一个实例”
# ＃制作一个KNeighborsClassifier对象的实例

knn = KNeighborsClassifier(n_neighbors=1)
type(knn)

# - 创建一个“知道”如何做K最近邻分类的对象，只是在等待数据
# - 对象的名称不重要
# - 在此步骤中可以指定调整参数（又名“超参数”）
# - 未指定的所有参数都设置为其默认值

print knn

#**步骤3：**将模型与数据（也称为“模型训练”）相匹配
# - 模型是在我们的“训练数据”中“学习”X和y之间的关系，
# - 学习发生的过程因型号而异
# - 发生在原地
knn.fit(X, y)

# 一旦模型适合数据，就称为“拟合模型”
#
# **步骤4：**预测新观测的响应
#  - 新的观察被称为“样本外”数据
#  - 使用模型训练过程中学到的信息

#- 返回一个NumPy数组，我们跟踪数字“意味着什么”
#- 可以一次预测多个观测值

X_new = [[3, 5, 4, 2], [5, 4, 3, 2]]
knn.predict(X_new)

###调整一个KNN模型

#实例化模型（使用值K = 5）

knn = KNeighborsClassifier(n_neighbors=5)

# 用数据拟合模型
knn.fit(X, y)

# 预测新观察的响应
knn.predict(X_new)

knn.predict_proba(X_new)

##比较KNN和其他算法模型
# ** KNN的优势：**
# - 简单理解和解释
# - 模特训练很快
# - 可以用于分类和回归

# ** KNN的缺点：**
#  - 必须存储所有的训练数据
#  当n很大时，预测阶段会变慢
# - 敏感的不相关的功能
#  - 对数据的规模敏感
#  - 准确性（一般）不与最好的监督学习方法竞争

第五课、KNN算法分析NBA球员数据

# # KNN运动与NBA球员数据

# **进球：**使用助攻，抢断，盖帽，失误和个人犯规来预测球员的位置
# 步骤1：将数据读入Pandas
# 将数据读入DataFrame

import pandas as pd

nba = pd.read_csv('NBA_players_2015.csv',index_col=0)

# 检查列
nba.columns

# 检查位置
nba.pos.value_counts()

# 步骤2：创建X和y
# ＃使用以下功能：助攻，抢断，盖帽，失误，个人犯规
#
# ＃将位置映射到数字
nba['pos_num'] = nba.pos.map({'C':0, 'F':1, 'G':2})

#创建特征矩阵（X）
feature_cols = ['ast', 'stl', 'blk', 'tov', 'pf']
X = nba[feature_cols]

# 替代方法来创建X
X = nba.loc[:, 'ast':'pf']

# 创建响应向量（y）
y = nba.pos_num

# ## 第三步：训练一个KNN模型（K = 5）

# import class
from sklearn.neighbors import KNeighborsClassifier

# 用K = 5来实例化
knn = KNeighborsClassifier(n_neighbors=5)

# 符合数据
knn.fit(X, y)

#步骤4：预测玩家位置并计算每个位置的预测概率
# 预测这些统计数据：1次助攻，1次抢断，0次封盖，1次失误，2次犯规

# 创建一个列表来表示一个玩家
player = [[1, 1, 0, 1, 2]]

# 做一个预测
knn.predict(player)

# 计算预测概率
knn.predict_proba(player)

##步骤5：使用K = 50重复步骤3和4

# 重复K = 50
knn = KNeighborsClassifier(n_neighbors=50)
knn.fit(X, y)
knn.predict(player)

# 计算预测概率
knn.predict_proba(player)  #测试数据X的返回概率估计。

##奖金：探索功能，以决定哪些是预测性的

# 允许重复出现在笔记本中
import matplotlib.pyplot as plt

# 增加默认图形和字体大小以便于查看
plt.rcParams['figure.figsize'] = (6, 4)
plt.rcParams['font.size'] = 14

# 按位置分组的助攻说明
nba.groupby('pos').ast.describe().unstack()

# 按位置分组的助攻盒图
nba.boxplot(column='ast', by='pos')

# 按位置分组的助攻直方图
nba.hist(column='ast', by='pos', sharex=True)

第六课、模型评估

#**解决方案：**创建一个程序**估计**模型可能对样本外数据执行的可能性以及如何在模型之间进行选择。
#**注意：**这些程序可以与**任何机器学习模式**一起使用，不仅仅是KNN。
##评估程序＃1：训练和测试整个数据集
# ＃1.在整个数据集**上训练模型。
# ＃2.在相同数据集**上测试模型，并通过比较**预测**响应值和**真实**响应值来评估我们做得如何。

import pandas as pd
nba = pd.read_csv('NBA_players_2015.csv',index_col=0)

#将位置映射到数字
nba['pos_num'] = nba.pos.map({'C': 0, 'F': 1, 'G': 2})

# 创建特征矩阵（X）
feature_cols = ['ast', 'stl', 'blk', 'tov', 'pf']
X = nba[feature_cols]

# 创建响应向量（y）
y = nba.pos_num

### KNN (K=50)

from sklearn.neighbors import KNeighborsClassifier

# 实例化模型
knn = KNeighborsClassifier(n_neighbors=50)

# 在整个数据集上训练模型
knn.fit(X, y)

# 预测X中观测值的响应值（“测试模型”）
knn.predict(X)

# 存储预测的响应值
y_pred_class = knn.predict(X)

# 要评估一个模型，我们还需要一个**评估指标：**
#  - 用于量化模型性能的数字计算
#  - 适当的指标取决于你的问题的**目标**
# 分类问题最常见的选择：
#  **分类准确度**：正确预测的百分比（“奖励函数”，因为越高越好）
#  **分类错误**：错误预测的百分比（“损失函数”，因为越低越好）
# 在这种情况下，我们将使用分类准确性。

#计算分类准确度

from sklearn import metrics

print metrics.accuracy_score(y, y_pred_class)

#这被称为**训练准确性**，因为我们正在用我们用来训练模型的相同数据来评估模型。

# ### KNN (K=1)

knn = KNeighborsClassifier(n_neighbors=1)
knn.fit(X, y)
y_pred_class = knn.predict(X)
print metrics.accuracy_score(y, y_pred_class)

#在相同的数据上进行培训和测试的问题
#- 目标是估计**样本外的数据的可能性能**
#- 但是，最大限度地提高训练准确性奖励**过于复杂的模型**不一定一概而论
#- 不必要的复杂模型** overfit **训练数据：
#- 使用样本内数据进行测试时效果会很好
#- 对样本外的数据可能做得不好
#- 学习数据中的“噪音”，而不是“信号”
#- 来自Quora：[过度拟合的直觉解释是什么？]
# （http://www.quora.com/What-is-an-intuitive-explanation-of-overfitting/answer/Jessica-Su）
# **因此，训练准确性不是对样本外精确度的一个很好的估计。**
# ！[1NN分类图]（images / iris_01nn_map.png）
###评估程序＃2：培训/测试拆分
# ＃1.将数据集分成两部分：一个**训练集**和一个**测试集**。
# ＃2.在**训练集**上训练模型。
# ＃3.在**测试集**上测试模型，并评估我们做得如何。

# 这完成了什么？
# 模型可以在**不同的数据**上进行训练和测试（我们将测试数据视为超出样本数据）。
#  - 测试集的响应值是已知的，因此可以评估**预测。
# 这被称为**测试准确性**，因为我们正在模型训练中未使用的独立“测试集”上评估模型。

# **测试精度是训练精度比样本外性能更好的估计。**

def min_max(nums):
    smallest = min(nums)
    largest = max(nums)
    return [smallest, largest]

min_and_max = min_max([1, 2, 3])
print min_and_max
print type(min_and_max)

the_min, the_max = min_max([1, 2, 3])
print the_min
print type(the_min)
print the_max
print type(the_max)

# ### 了解`train_test_split`函数

from sklearn.cross_validation import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y)

# 在分裂之前
print X.shape

# 分裂后
print X_train.shape
print X_test.shape

#在分裂之前
print y.shape

# 分裂后
print y_train.shape
print y_test.shape

# ### 了解`random_state`参数

# 没有一个random_state参数
X_train, X_test, y_train, y_test = train_test_split(X, y)

#打印每个对象的第一个元素
print X_train.head(1)
print X_test.head(1)
print y_train.head(1)
print y_test.head(1)

#用一个random_state参数
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=99)

# 打印每个对象的第一个元素
print X_train.head(1)
print X_test.head(1)
print y_train.head(1)
print y_test.head(1)

# ###使用train/test 拆分程序（K = 1）

# 步骤1：将X和y分成训练和测试集（使用random_state进行再现）
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=99)

# 步骤2：在训练集上训练模型（使用K = 1）
knn = KNeighborsClassifier(n_neighbors=1)
knn.fit(X_train, y_train)

# 步骤3：在测试装置上测试模型，并检查准确性
y_pred_class = knn.predict(X_test)
print metrics.accuracy_score(y_test, y_pred_class)

# ### 重复K = 50

knn = KNeighborsClassifier(n_neighbors=50)
knn.fit(X_train, y_train)
y_pred_class = knn.predict(X_test)
print metrics.accuracy_score(y_test, y_pred_class)

# ！[偏差方差折衷]（images / bias_variance.png）
####比较测试精度和零精度
# 空值准确度是通过**始终预测最频繁的班级**所能达到的准确度。 这是您可能想要测量您的分类模型的基准。
#检查class分布
y_test.value_counts()

#计算null精度
y_test.value_counts().head(1) / len(y_test)

# ### 寻找K的“最佳”值

# 计算K = 1到100的TRAINING ERROR和TESTING ERROR

k_range = range(1, 101)
training_error = []
testing_error = []

for k in k_range:
    # 用当前K值实例化模型
    knn = KNeighborsClassifier(n_neighbors=k)

    # 计算训练错误
    knn.fit(X, y)
    y_pred_class = knn.predict(X)
    training_accuracy = metrics.accuracy_score(y, y_pred_class)
    training_error.append(1 - training_accuracy)

    # 计算测试错误
    knn.fit(X_train, y_train)
    y_pred_class = knn.predict(X_test)
    testing_accuracy = metrics.accuracy_score(y_test, y_pred_class)
    testing_error.append(1 - testing_accuracy)

# 让情节出现在笔记本中
import matplotlib.pyplot as plt

plt.style.use('fivethirtyeight')

# 创建K的DataFrame，训练错误和测试错误
column_dict = {'K': k_range, 'training error': training_error, 'testing error': testing_error}
df = pd.DataFrame(column_dict).set_index('K').sort_index(ascending=False)
df.head()

#绘制K（HIGH TO LOW）和TESTING ERROR之间的关系
df.plot(y='testing error')
plt.xlabel('Value of K for KNN')
plt.ylabel('Error (lower is better)')

# 找到最小的测试误差和相关的K值
df.sort_values('testing error').head()

# 替代方法
min(zip(testing_error, k_range))

#我们可以得出什么结论？
#在使用这些特征的数据集上使用KNN时，K **的**最佳值可能在14左右。
#给定**未知玩家**的统计数据，我们估计我们能够在74％的时间内正确地预测他的位置。
###训练错误与测试错误
#绘制K（高到低）和TRAINING ERROR和TESTING ERROR之间的关系
df.plot()
plt.xlabel('Value of K for KNN')
plt.ylabel('Error (lower is better)')

# **培训错误**随着模型复杂度的增加而降低（K值越低）
# **测试错误**在最佳模型复杂度下最小化
# [偏差方差折衷]（images / training_testing_error.png）
###对超出样本的数据进行预测
#根据一个（真正）未知的球员的统计数据，我们如何预测他的位置？
#用已知最好的参数实例化模型

knn = KNeighborsClassifier(n_neighbors=14)

# 为什么重新训练模型用X和Y（不是X_train和y_train）？
knn.fit(X, y)

# 对样本外观察进行预测
knn.predict([[1, 1, 0, 1, 2]])

###traning/test拆分的缺点？
#如果`train_test_split`函数以不同的方式分割数据会发生什么？ 我们会得到和以前一样的确切结果吗？
#为random_state尝试不同的值

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=98)
knn = KNeighborsClassifier(n_neighbors=50)
knn.fit(X_train, y_train)
y_pred_class = knn.predict(X_test) 
print metrics.accuracy_score(y_test, y_pred_class)

#测试准确性是样本外精度的**高方差估计**
# ** K倍交叉验证**克服了这个限制，并提供了更可靠的估计
# - 但是，train/test 拆分仍然是有用的，因为它**的灵活性和速度*

教程来源：https://github.com/justmarkham/DAT8 作者：Kevin Markham

#数据链接：https://pan.baidu.com/s/1o8eLDEi 密码：i1zx

# 本实战教程是python2.7版本

# 涉及的API
# http://pandas.pydata.org/pandas-docs/stable/genindex.html

# https://docs.python.org/2/library/index.html
#############
# "第一课：操作TSV文件"
############

'''
第1部分：用csv.reader（）读入文件，并将其存储在名为“file_nested_list”的对象中。
提示：这是一个TSV文件，需要告诉csv.resader（）如何处理它。

'''

import csv

# 指定分隔符是制表符，TSV为用制表符tab分隔的文本文件。 CSV为用逗号,分隔的文本文件。
with open('chipotle.tsv', mode='rU') as f:
    file_nested_list = [row for row in csv.reader(f, delimiter='\t')]

'''
第2部分：将'file_nested_list'分隔为'header'和'data'。
'''

header = file_nested_list[0]
data = file_nested_list[1:]   #从第二个元素开始截取列表，

'''
第3部分：计算订单的平均价格。
提示：检查数据以查看“数量”列是否与此计算相关。
提示：仔细想想最简单的方法来做到这一点！
'''

#统计唯一的order_id的数量
#注意：你可以认为这是1834年，因为这是最大的order_id，但最好检查一下
num_orders = len(set([row[0] for row in data]))     # 1834
data

# 对data列表遍历，因为里面还有列表，
# 因为要剥离美元符号和尾随空间，所以读data里面的列表取索引4的值，然后在对里面从索引1处到结尾的那个数字，进行取值，
#创建一个价格清单
prices = [float(row[4][1:-1]) for row in data]

prices

# calculate the average price of an order and round to 2 digits
round(sum(prices) / num_orders, 2)      # $18.81    round对浮点数进行近似取值，保留几位小数

'''
第4部分：创建一个销售的所有独特的苏打水和软饮料的清单（或集合）。
注意：只要看“'Canned Soda' 和'Canned Soft Drink'”，不要忽略其他饮料，如“Izze”。
'''

#如果“item_name”包含“Canned”，则将“choice_description”附加到“sodas”列表中
sodas = []
for row in data:
    if 'Canned' in row[2]:   #因为item_name在data列表中索引是2
        sodas.append(row[3][1:-1])      # 去掉括号，choice_description在data列表中是索引3，去长度索引为1，到结尾的字符。

# 等价上面，一个效果！
sodas = [row[3][1:-1] for row in data if 'Canned' in row[2]]

# 创建一个无重复的列表。
unique_sodas = set(sodas)
unique_sodas

'''
第5部分：计算每个卷饼的平均配料数量。
注意：我们忽略“数量”列来简化此任务。
'''

# 总的卷饼和配料
burrito_count = 0    #卷饼数量
topping_count = 0    #配料数量

# 通过计算逗号并加1来计算配料的数量
# note: x += 1 is equivalent to x = x + 1
for row in data:
    if 'Burrito' in row[2]:
        burrito_count += 1
        topping_count += (row[3].count(',') + 1)

# calculate the average topping count and round to 2 digits
round(topping_count / float(burrito_count), 2)      # 5.40

'''
第6部分：创建一个字典，其中的键代表芯片订单和该值表示订单总数。
预计产量：{'芯片和烤辣椒玉米莎莎'：18，...}
注意：请考虑“数量”栏！
可选：了解如何使用“defaultdict”来简化代码。
'''

#创建一个空字典
chips = {}

# 如果芯片顺序不在字典中，则添加一个新的键/值对
# 如果芯片顺序已经在字典中，则更新该键的值
for row in data:
    if 'Chips' in row[2]:
        if row[2] not in chips:
            chips[row[2]] = int(row[1])     # 这是一个新键，所以创建键/值对，row[1]代表订单（quantity）数量，
        else:
            chips[row[2]] += int(row[1])    # 这是一个现有的键值，所以当遍历到新的chips没有在chips字典里时候，所以增加值。

# defaultdict 函数可以帮你检查key时候已经存在的麻烦~，跟上面的for 代码一样！简化了一层if！
from collections import defaultdict
dchips = defaultdict(int)
for row in data:
    if 'Chips' in row[2]:
        dchips[row[2]] += int(row[1])

chips

#############
# "第二课：pandas统计操作csv文件"
############

'''
基本水平
'''

import pandas as pd
import matplotlib.pyplot as plt

# 读入“imdb_1000.csv”并将其存储在名为movies的DataFrame中
movies = pd.read_csv('imdb_1000.csv')

# 检查行数和列数
movies.shape

# 检查每列的数据类型
movies.dtypes

# 计算平均电影持续时间，mean是均值统计函数！  计算出这1千部电影，电影平均持续时间大约为120分钟！
movies.duration.mean()
# 按持续时间排序DataFrame以查找最短和最长的电影
movies.sort_values('duration').head(1)   #时间最短的电影
movies.sort_values('duration').tail(1) #时间最长的电影：Kenneth Branagh'

# 创建一个持续时间的直方图，选择一个“适当”的箱子数量
movies.duration.plot(kind='hist', bins=20)

# 使用箱形图来显示相同的数据
movies.duration.plot(kind='box')

'''
中级水平
'''

# 统计每个内容分级的电影数量,分级字(content_rating)
movies.content_rating.value_counts()  #返回包含唯一值计数的对象,默认以降序排列

# 使用可视化显示相同的数据，包括标题和x和y标签
plt.xlabel('Content Rating')
plt.ylabel('Number of Movies')
movies.content_rating.value_counts().plot(kind='bar', title='Top 1000 Movies by Content Rating')

# 将 NOT RATED, APPROVED, PASSED, GP 评级 都转换为"UNRATED"，初始状态为38个！
movies.content_rating.replace(['NOT RATED', 'APPROVED', 'PASSED', 'GP'], 'UNRATED', inplace=True)

movies.content_rating.value_counts() # 验证

# 将 'X', 'TV-MA评级 都转换为"NC-17"!
movies.content_rating.replace(['X', 'TV-MA'], 'NC-17', inplace=True)

#计算每列中缺失值的数量
movies.isnull().sum()

# 如果缺少值：检查它们，然后用“合理的”值填充它们
movies[movies.content_rating.isnull()]
movies.content_rating.fillna('UNRATED', inplace=True)

# 计算电影2小时以上的平均星级评分，
# 并与短于2小时的电影的平均星级进行比较
movies[movies.duration >= 120].star_rating.mean()
movies[movies.duration < 120].star_rating.mean()

# 使用可视化来检测持续时间和星级之间是否存在关系
movies.plot(kind='scatter', x='duration', y='star_rating', alpha=0.2)

# 计算每个类型的平均持续时间
movies.groupby('genre').duration.mean()

'''
先进水平
'''

# 可视化内容评级和持续时间之间的关系
movies.boxplot(column='duration', by='content_rating')
movies.hist(column='duration', by='content_rating', sharex=True)

# 确定最高评级的电影（按星级评分）为每个流派
movies.sort_values('star_rating', ascending=False).groupby('genre').title.first()
movies.groupby('genre').title.first()   # 因为DataFrame已经按照星级评分

# 检查是否有多个具有相同标题的电影，如果是，确定它们是否实际上是重复的
dupe_titles = movies[movies.title.duplicated()].title
movies[movies.title.isin(dupe_titles)]  # 将重复的标题显示出来

# 计算每个流派的平均星级，但只包括至少有10部电影的流派

# 选项1：手动创建相关流派的列表，然后使用该列表进行过滤
movies.genre.value_counts()
top_genres = ['Drama', 'Comedy', 'Action', 'Crime', 'Biography', 'Adventure', 'Animation', 'Horror', 'Mystery']
movies[movies.genre.isin(top_genres)].groupby('genre').star_rating.mean()

# 选项2：通过保存value_counts然后过滤，自动创建相关流派的列表
genre_counts = movies.genre.value_counts()
top_genres = genre_counts[genre_counts >= 10].index

movies[movies.genre.isin(top_genres)].groupby('genre').star_rating.mean()

# 选项3：计算所有流派的平均星级评分，然后使用布尔序列进行筛选
movies.groupby('genre').star_rating.mean()[movies.genre.value_counts() >= 10]

# 选项4：通过计数和平均值聚合，然后使用计数进行过滤
genre_ratings = movies.groupby('genre').star_rating.agg(['count', 'mean'])
genre_ratings[genre_ratings['count'] >= 10]

#############
# "第三课：鸢尾花数据用来预测种类
############

# 增加默认图形和字体大小以便于查看
plt.rcParams['figure.figsize'] = (8, 6)
plt.rcParams['font.size'] = 14

# ## Task 1
#
#将鸢尾花数据读入熊猫数据框，包括列名。

#定义列名称列表（作为字符串）(sepal)萼片,(petal)花瓣,(species)种类
col_names = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'species']

# 定义从中检索数据的URL（以字符串形式）
url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data';;

# 检索CSV文件并添加列名称
iris = pd.read_csv(url, header=None, names=col_names)

# ## Task 2
#
# 收集有关数据的一些基本信息。

iris.shape

iris.head()

iris.dtypes

iris.describe()

iris.species.value_counts()

iris.isnull().sum()

# ## Task 3
#
# 使用排序，拆分应用组合和/或可视化来查找物种之间的差异。.

# ### sorting

# 通过petal_width对DataFrame进行排序并显示NumPy数组

iris.sort_values('petal_width').values

# ### 拆分申请，结合
# 是按物种分组，萼片的长度的统计均值。
iris.groupby('species').sepal_length.mean()

# 所有数字列的均值按物种分组
iris.groupby('species').mean()

#描述按物种分组的所有数字列
iris.groupby('species').describe()  #生成描述性统计数据，汇总数据集分布的中心趋势。

# ### 可视化

# 花瓣宽度的直方图按物种分组
iris.hist(column='petal_width', by='species', sharex=True)

# 根据物种分组的花瓣宽度箱形图
iris.boxplot(column='petal_width', by='species')

# 按种类分组的所有数字列的框图
iris.boxplot(by='species')

# 将物种映射到数值，以便可以按物种染色
iris['species_num'] = iris.species.map({'Iris-setosa': 0, 'Iris-versicolor': 1, 'Iris-virginica': 2})

# 替代方法
iris['species_num'] = iris.species.factorize()[0]

# petal_length（花瓣长度） vs petal_width（花瓣宽度）的散点图
iris.plot(kind='scatter', x='petal_length', y='petal_width', c='species_num', colormap='brg')

# 所有特征的物种的散布矩阵
pd.scatter_matrix(iris.drop('species_num', axis=1), c=iris.species_num, figsize=(12, 10))

# ## Task 4
#
# 写下一套可用于根据鸢尾花数据测量结果预测物种的规则。.

# 定义一个代表花瓣区域的新特征（“特征工程”）
iris['petal_area'] = iris.petal_length * iris.petal_width

# 描述花瓣区域按物种分组
iris.groupby('species').petal_area.describe().unstack()#旋转（必须是分层的）索引标签的级别，返回具有新标签级别的DataFrame

# 根据物种分组的花瓣区域箱形图
iris.boxplot(column='petal_area', by='species')

    # 只显示7至9的花瓣区域的鸢尾花
iris[(iris.petal_area > 7) & (iris.petal_area < 9)].sort_values('petal_area')

# 我的一套预测物种的规则：
#
# - 如果petal_area小于2，则预测** setosa **
# - - 如果petal_area小于7.4，则预测**versicolor*.
# -否则，预测** virginica ***.

# ## Bonus
#
# 定义一个接受一行数据并返回预测物种的函数。 然后，使用该函数对所有现有数据行进行预测，并检查预测的准确性。
#给出一行数据，返回预测的species_num（0/1/2）
def classify_iris(row):
    # c计算petal_area
    petal_area = row[2] * row[3]

    # 根据上述规则预测物种
    if petal_area < 2:
        prediction = 'setosa'
    elif petal_area < 7.4:
        prediction = 'versicolor'
    else:
        prediction = 'virginica'

    # 将物种名称映射到数值
    species_to_num = {'setosa': 0, 'versicolor': 1, 'virginica': 2}

    # 返回该值
    return species_to_num[prediction]

# 打印第一行
iris.iloc[0, :]

# 打印最后一行
iris.iloc[149, :]

# 测试第一行和最后一行的功能
print classify_iris(iris.iloc[0, :])
print classify_iris(iris.iloc[149, :])

# 对所有行进行预测并将其存储在DataFrame中
iris['prediction'] = [classify_iris(row) for index, row in iris.iterrows()]

# 计算正确预测的百分比
sum(iris.species_num == iris.prediction) / 150.

##教程来源：https://github.com/justmarkham/DAT8   作者：Kevin Markham

# -*-coding:utf-8-*-

# ---------------------
# Chapter 1 - Reading from a CSV.ipynb
# ---------------------

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

"""
    1.1 Reading data from a csv file
"""
broken_df = pd.read_csv('data/bikes.csv')

# seq为分隔符，encoding为编码，index_col为索引列编号，dayfirst为日期格式，parse_dates为日期解析
fixed_df = pd.read_csv('data/bikes.csv', sep=';', encoding='latin1', parse_dates=['Date'], dayfirst=True,
                       index_col='Date')
fixed_df[:3]
type(fixed_df)  # <class 'pandas.core.frame.DataFrame'>

"""
    1.2 Selecting a column
"""
fixed_df['Berri 1']

"""
    1.3 Plotting a column
"""
fixed_df['Berri 1'].plot()
plt.show()

fixed_df.plot(figsize=(15, 10))
plt.show()

"""
    1.4 Putting all that together
"""
df = pd.read_csv('data/bikes.csv', sep=';', encoding='latin1', parse_dates=['Date'], dayfirst=True, index_col='Date')
df['Berri 1'].plot()
plt.show()

# ---------------------
# Chapter 2 - Selecting data & finding the most common complaint type.ipynb
# ---------------------

pd.set_option('display.line_width', 5000)
pd.set_option('display.max_columns', 60)
complaints = pd.read_csv('data/311-service-requests.csv')
complaints.head()
complaints['Complaint Type']
complaints[:3]
complaints['Complaint Type'][:3]
complaints[:3]['Complaint Type']
complaints[['Complaint Type', 'Borough']]
complaints[['Complaint Type', 'Borough']][:10]
complaints_counts = complaints['Complaint Type'].value_counts()  # 计算各个元素的数量
complaints_counts[:10]
complaints_counts[:10].plot(kind='bar')
plt.show()

# ---------------------
# Chapter 3 - Which borough has the most noise complaints (or, more selecting data).ipynb
# ---------------------

complaints = pd.read_csv('data/311-service-requests.csv')  # 读取csv文件
complaints.head()
complaints[:5]

'''
    3.1 Selecting only noise complaints
'''
noise_complaints = complaints[complaints['Complaint Type'] == 'Noise - Street/Sidewalk']
noise_complaints[:3]
complaints['Complaint Type'] == 'Noise - Street/Sidewalk'  # 返回True False

is_noise = complaints['Complaint Type'] == 'Noise - Street/Sidewalk'
in_brooklyn = complaints['Borough'] == 'BROOKLYN'
complaints[is_noise & in_brooklyn][:5]

complaints[is_noise & in_brooklyn][['Complaint Type', 'Borough', 'Created Date', 'Descriptor']][:10]
'''
    3.2 A digression about numpy arrays
'''
pf = pd.Series([1, 2, 3])
pf
pf.values
pf.index
nf = np.array([1, 2, 3])
nf
nf != 2
nf[nf != 2]

'''
    3.3 So, which borough has the most noise complaints?
'''
is_noise = complaints['Complaint Type'] == "Noise - Street/Sidewalk"
noise_complaints = complaints[is_noise]
noise_complaints['Borough'].value_counts()

noise_complaint_counts = noise_complaints['Borough'].value_counts()
complaint_counts = complaints['Borough'].value_counts()
noise_complaint_counts / complaint_counts
noise_complaint_counts / complaint_counts.astype(float)
(noise_complaint_counts / complaint_counts.astype(float)).plot(kind='bar')
plt.show()

# ---------------------
# Chapter 4: Find out on which weekday people bike the most with groupby and aggregate
# ---------------------

bikes = pd.read_csv('data/bikes.csv', sep=';', encoding='latin1', index_col='Date', parse_dates=['Date'],
                    dayfirst=True)
bikes.head()

bikes['Berri 1'].plot()  # 绘制曲线
plt.show()

berri_bikes = bikes[['Berri 1']].copy()  # 将某一列的数据复制出来，单独为一列
berri_bikes[:5]
berri_bikes.index
berri_bikes.index.day
berri_bikes.index.weekday
berri_bikes.loc[:, 'weekday'] = berri_bikes.index.weekday
berri_bikes[:5]

"""
    4.2 Adding up the cyclists by weekday
"""
"""
    使用DataFrames中的.groupby()方法进行分组，并计算每一组的数量和
"""
weekday_counts = berri_bikes.groupby('weekday').aggregate(sum)
weekday_counts

weekday_counts.index = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday']
weekday_counts

weekday_counts.plot(kind='bar')
plt.show()

"""
    4.3 Putting it together
"""
"""
    所有代码汇总
"""
bikes = pd.read_csv('data/bikes.csv', sep=';', encoding='latin1', index_col='Date', dayfirst=True,
                    parse_dates=['Date'])
berri_bikes = bikes[['Berri 1']].copy()
berri_bikes.loc[:, 'weekday'] = berri_bikes.index.weekday

weekday_counts = berri_bikes.groupby('weekday').aggregate(sum)
weekday_counts.index = ['Monday', 'Tuesday', 'Wednesday', 'Thursday', 'Friday', 'Saturday', 'Sunday']
weekday_counts.plot(kind='bar')
plt.show()

"""
    分析：
        主要是计算时间，分组处理一周时间，将每周对应的数量加到对应的天上
    方法：
        1、csv数据的读取
        2、列数据的复制
        3、将数据按照一周来进行划分
        4、按照一周进行分组处理数据，修改索引
        5、直方图展示
"""

# ---------------------
# Chapter 5: Combining dataframes and scraping Canadian weather data
# ---------------------

'''
    Summary
'''
weather_2012_final = pd.read_csv('data/weather_2012.csv', index_col='Date/Time')
weather_2012_final.head()
weather_2012_final['Temp (C)'].plot(figsize=(15, 6))
plt.show()

'''
    5.1 Downloading one month of weather data
'''

url_template = "http://climate.weather.gc.ca/climate_data/bulk_data_e.html?format=csv&stationID=31407&Year={year}&Month={month}&timeframe=1&submit=%E4%B8%8B%E8%BD%BD%E6%95%B0%E6%8D%AE";;
url = url_template.format(month=3, year=2012)
weather_mar2012 = pd.read_csv(url, skiprows=15, index_col='Date/Time', parse_dates=True, encoding='latin1', header=0)

weather_mar2012

weather_mar2012['Temp (°C)'].plot(figsize=(15, 5))  # 图形展示温度变化情况
plt.show()

weather_mar2012.columns = [
    u'Year', u'Month', u'Day', u'Time', u'Data Quality', u'Temp (°C)',
    u'Temp Flag', u'Dew Point Temp (C)', u'Dew Point Temp Flag',
    u'Rel Hum (%)', u'Rel Hum Flag', u'Wind Dir (10s deg)', u'Wind Dir Flag',
    u'Wind Spd (km/h)', u'Wind Spd Flag', u'Visibility (km)', u'Visibility Flag',
    u'Stn Press (kPa)', u'Stn Press Flag', u'Hmdx', u'Hmdx Flag', u'Wind Chill',
    u'Wind Chill Flag', u'Weather']

# weather_mar2012 = weather_mar2012.dropna(axis=1, how='any')  # drop the column if any value is null 删除空列
weather_mar2012[:5]

# weather_mar2012 = weather_mar2012.drop(['Year', 'Month', 'Day', 'Time'], axis=1)

weather_mar2012

'''
    5.2 Plotting the temperature by hour of day
'''
temperatures = weather_mar2012[['Temp (°C)']].copy()
temperatures.head()

temperatures.loc[:, 'Hour'] = weather_mar2012.index.hour
temperatures.groupby('Hour').aggregate(np.median).plot()
plt.show()

'''
    5.3 Getting the whole year of data
'''

def download_weather_month(year, month):
    if month == 1:
        year += 1
    url = url_template.format(year=year, month=month)
    weather_data = pd.read_csv(url, skiprows=15, index_col='Date/Time', parse_dates=True, header=True)
    weather_data = weather_data.dropna(axis=1)
    weather_data.columns = [col.replace('\xb0', '') for col in weather_data.columns]
    weather_data = weather_data.drop(['Year', 'Day', 'Month', 'Time', 'Data Quality'], axis=1)
    return weather_data

download_weather_month(2012, 1)[:5]

data_by_month = [download_weather_month(2012, i) for i in range(1, 12)]  # 所有月份

weather_2012 = pd.concat(data_by_month)
weather_2012

'''
    5.4 Saving to a CSV
'''
weather_2012.to_csv('data/weather_2012.csv')

# ---------------------
# Chapter 6 - String Operations- Which month was the snowiest.ipynb
# ---------------------

weather_2012 = pd.read_csv('data/weather_2012.csv', encoding='latin1', index_col='Date/Time', parse_dates=True)
weather_2012.head()

'''
    6.1 String operations
'''
weather_description = weather_2012['Weather']  # 取出天气那一列

# 天气是否为snow，若天气为snow，返回True,否则返回False；str将类型转换成字符串形式，利于字符串的匹配、替换和截取
is_snowing = weather_description.str.contains('Snow')

is_snowing = is_snowing.astype(float)

is_snowing.plot()  # 将一年中下雪的天全部显示出来
plt.show()

'''
    6.2 Use resampling to find the snowiest month 寻找下雪最多的月份
'''
'''
    每个月份的平均气温，可以使用resample()方法来实现
'''
weather_2012['Temp (C)'].resample('M', how=np.median).plot(kind='bar')  # 平均气温
plt.show()

# 将天气情况使用0和1来表示，若为snow，也就是True，则返回1，否则用0表示
is_snowing.astype(float)[:10]  # astype用于类型转换，bool类型转换为float类型

# 使用resample()查找出每个月下雪的可能性，用百分比来表示
is_snowing.astype(float).resample('M', how=np.mean)  # 'M':表示按月的时间频率

is_snowing.astype(float).resample('M', how=np.mean).plot(kind='bar')
plt.show()

'''
    6.3 Plotting temperature and snowiness stats together  温度和下雪一起分析
'''
temperature = weather_2012['Temp (C)'].resample('M', how=np.median)  # 平均温度
is_snowing = weather_2012['Weather'].str.contains('Snow')
snowiness = is_snowing.astype(float).resample('M', how=np.mean)  # 下雪比例

temperature.name = 'Temperature'
snowiness.name = 'Snowiness'

stats = pd.concat([temperature, snowiness], axis=1)
stats
stats.plot(kind='bar')  # 下雪比例在图示中显示太小，所以这里不合理
plt.show()

stats.plot(kind='bar', subplots=True, figsize=(15, 10))  # 将两张图放到一个平面上，分开放，这样就能合理的进行图像展示
plt.show()

# subplot()作用：将多个figure放到一个平面上

# ---------------------
# Chapter 7 - Cleaning up messy data.ipynb  清理垃圾数据
# ---------------------

requests = pd.read_csv('data/311-service-requests.csv')
# requests.head()

'''
    7.1 How do we know if it's messy?
'''
zip = requests['Incident Zip'].unique()  # unique()用于查看所有的值
# zip
'''
    zip中存在的问题：
        1、数据类型问题，有些是字符串型，有些是浮点型
        2、有一些值不存在nan
        3、有些值不正确 83  29616-0759
        4、有N/A值，pandas不能够识别，'N/A','NO CLUE'
    处理方法：
        1、使'N/A','NO CLUE'变成规则的nan
        2、使所有格式都变成字符串
'''

'''
    7.3 Fixing the nan values and string/float confusion
'''
na_value = ['N/A', 'NO CLUE', 'O', 'nan']
requests = pd.read_csv('data/311-service-requests.csv', na_values=na_value, dtype={'Incident Zip': str})
# 读取csv文件时，将异常值设置为空值，将数据类型全部转换为字符串类型
zip = requests['Incident Zip'].copy()
# zip.unique()

'''
    7.4 What's up with the dashes? 处理数字之间的横杠29616-0759
'''
row_with_dashs = requests['Incident Zip'].str.contains('-').fillna(False)  # 将带横杠的全部提取出来
# len(requests[row_with_dashs])
# requests[row_with_dashs]

requests['Incident Zip'][row_with_dashs] = np.nan  # 将带横杠的全部转换为空值
# requests['Incident Zip'].unique()

long_zip_codes = requests['Incident Zip'].str.len() > 5
# requests['Incident Zip'][long_zip_codes].unique()

requests['Incident Zip'] = requests['Incident Zip'].str.slice(0, 5)  # slice()获取字符串的指定长度
# requests['Incident Zip'] = requests['Incident Zip'].str[0:5]
# requests['Incident Zip'].unique()

# requests[requests['Incident Zip']] == '00000'

zero_zips = requests['Incident Zip'] == '00000'
requests.loc[zero_zips, 'Incident Zip'] = np.nan

unique_zips = requests['Incident Zip'].unique()
unique_zips.sort()  # 排序
unique_zips

zips = requests['Incident Zip']
is_close = zips.str.startswith('0') | zips.str.startswith('1')  # zip以0或1开头
is_far = ~(is_close) & zips.notnull()

zips[is_far]

requests[is_far][['Incident Zip', 'Descriptor', 'City']].sort('Incident Zip')

requests['City'].str.upper().value_counts()  # 城市名转换为大写的，并且统计城市的数量

'''
    7.5 Putting it together
'''
# 异常值处理及csv文件的读取
na_values = ['NO CLUE', 'N/A', '0']
requests = pd.read_csv('data/311-service-requests.csv',
                       na_values=na_values,
                       dtype={'Incident Zip': str})

# 将邮政编码的位数固定为5位
def fix_zip_codes(zips):
    zips = zips.str.slice(0, 5)

    zero_zips = zips == '00000'
    zips[zero_zips] = np.nan

    return zips

requests['Incident Zip'] = fix_zip_codes(requests['Incident Zip'])
requests['Incident Zip'].unique()

# ---------------------
# Chapter 8 - How to deal with timestamps.ipynb
# ---------------------

'''
    8.1 Parsing Unix timestamps
'''
popcon = pd.read_csv('data/popularity-contest', sep=' ')
#  popcon.head()
popcon.columns = ['atime', 'ctime', 'package-name', 'mru-program', 'tag']
#  popcon[:5]
popcon['atime'].dtype

popcon['atime'] = popcon['atime'].astype(int)
#  popcon['atime'][:5]
# popcon['ctime'] = popcon['ctime'].astype(int)
popcon['atime'] = pd.to_datetime(popcon['atime'])
# popcon['ctime'] = pd.to_datetime(popcon['ctime'], unit='s')
#  popcon['atime'][:5]

popcon = popcon[popcon['atime'] > '1970-01-01']
nonlibraries = popcon[~popcon['package-name'].str.contains('lib')]
nonlibraries.sort('ctime', ascending=False)[:10]

# ---------------------
# Chapter 9 - Loading data from SQL databases.ipynb
# ---------------------

import sqlite3
import pandas as pd
"""
    pandas can read from HTML,JSON,SQL,EXCEL,HDF5,Stata, and a few other things.

    Read data from a SQL database using the pd.read_sql function.

    read_sql take 2 arguments: a SELECT statement, and s database connection object.

    This is great because it means you can read from any kind of SQL database,
    it doesn't matter if it's MySQL,SQLite,PostgreSQL,or something else.
"""
"""
    9.1 Reading data from SQL databases  读取数据
"""
con = sqlite3.connect('../data/weather_2012.sqlite')
df = pd.read_sql('select * from weather_2012 LIMIT 3', con, index_col='id')  # 设置id索引
#  df
df = pd.read_sql('select * from weather_2012 LIMIT 3', con, index_col=['id', 'date_time'])  # 设置双重索引
#  df

"""
    9.2 Writing to a SQLite database  写入数据
"""
# weather_df = pd.read_csv('../data/weather_2012.csv')
# con = sqlite3.connect('../data/test_db.sqlite')
# con.execute('drop table if exists weather_2012')
# weather_df.to_sql('weather_2012', con)

con = sqlite3.connect('../data/test_db.sqlite')
df = pd.read_sql('select * from weather_2012 LIMIT 3', con, index_col='index')
#  df

con = sqlite3.connect('../data/test_db.sqlite')
df = pd.read_sql('select * from weather_2012 order by Weather LIMIT 3', con)
df

"""
    sqlite3 database:连接数据库-->sqlite3.connect()
    PostgreSQL database:连接数据库-->psycopg2.connect()
    MySQL database:连接数据库-->MySQLdb.connect()
"""

"""
    9.3 Connecting to other kinds of database
"""
import MySQLdb
con = MySQLdb.connect(host='localhost', db='test')

import psycopg2
con = psycopg2.connect(host='localhost')