【机器学习(一)】预测建模

1.准备预测模型

◾ 机器学习的特征
◾ 什么是模型?
◾ 如何对模型进行评估?
◾ 工具:
第一个模型:k近邻(k-Nearest Neighbor)模型
第一个R包:caret

什么是“机器学习”(machine learning)?

机器学习:掌握一些算法,通过训练数据输入获得模型,利用模型学习新的数据,获得对预测性能的改善。
机器学习

1.1模型

模型:定量地或定性地描述系统各变量之间的相互关系或因果关系的一系列数学方程。

模型获得的特征一
◾ 根据数据来构建和选择模型
◾ 不以数学推理或逻辑归纳的手段从已知事实获得模型

例如:预测降雨量,收集有限地点的历史降雨量数据(降雨量,及高度、经度、纬度),来预测没有收集或者无法收集任何数据的地点的降雨量。

1.1.1 从数据中学习

模型获得的特征二
◾ 根据模型数据创建模型来描述现象的过程中,必然存在某些随机性的来源

模型包含:随机成分+确定成分

Iris Data Set(鸢尾属植物数据集)

数据集合包含三种鸢尾属植物[setosa,versicolor, virginica]。
每一种植物有50株植物,包含长度和宽度的数据(cm),共150株植物的数据[sepal:萼片;petal:花瓣]。

head(iris)
iris

令𝐷 = { 𝑥1 , 𝑥2, ⋯ , 𝑥𝑚} 表示包含 𝑚 个示例的数据集,每个示例由 𝑑 个属性描述, 则每个示例 𝑥𝑖 = (𝑥𝑖1, 𝑥𝑖2, ⋯ , 𝑥𝑖𝑑) 是 d维样本空间𝒳中的一个向量,𝑥𝑖 ∈ 𝒳,其中𝑥𝑖𝑗是示例𝑥𝑖在第𝑗个属性上的取值,𝑑 称为样本 𝑥𝑖的“维数”(dimensionality)。


library(ggplot2)
ggplot(iris, aes(Petal.Length)) + geom_histogram(bins = 30)
# bins指定分组数目(格子),默认值是30,如果未设置,将显示警告消息。
ggplot(iris, aes(Petal.Length, Petal.Width, color=Species)) + 
geom_point() + theme(legend.position = "bottom") 
# 图例显示在底部
预测建模的过程

创建一个函数来预测数量值,并使它产生的相对目标函数的误差最小化。

为什么我们总是很难找到目标函数?
◾ 可简化误差(reducible error):实际模型是非线性的用线性模型
◾ 不可简化误差(irreducible error):噪声,或者由于数据集合太小

1.1.2 模型的核心组成部分

● 带有需要调优参数的一组方程:线性回归、神经网络、支持向量机等有特定的参数化方程
● 代表我们建模所针对系统或过程的一些数据
● 描述该模型拟合优度的一个概念
● 更新参数以改善模型拟合优度的一种方法

kNN模型

K-nearest Neighbor
如果一个样本在特征空间中的K个最相似(即特征空间中最邻近)的样本中的大多数属于某一个类别,则该样本也属于这个类别。该方法在定类决策上只依据最邻近的一个或者几个样本的类别来决定待分样本所属的类别。


◾ 给定𝑋后先估计𝑌的条件分布,然后将一个给定的观测分类到估计分布概率最大的类别中

KNN与K-means的区别:
KNN属于监督学习,类别是已知的,通过对已知分类的数据进行训练和学习,找到这些不同类的特征,再对未分类的数据进行分类。
K-means属于非监督学习,事先不知道数据会分为几类,通过聚类分析将数据聚合成几个群体。聚类不需要对数据进行训练和学习。

用kNN模型来预测该花属于哪个品种?

第一步:确定新样本的k个最邻近样本

第二步:找到新样本到iris数据集里每个点的距离,然后对结果排序

apply()函数只能用于处理矩阵类型的数据。
apply()一般有三个参数:
第一个参数代表矩阵对象;
第二个参数代表要操作矩阵的维度,1表示对行进行处理,2表示对列进行处理;
第三个参数是处理数据的函数。apply会分别一行或一列处理该矩阵的数据。

# 新建一个样本
new_sample <- c(4.8, 2.9, 3.7, 1.7)
names(new_sample)<-c("Sepal.Length","Sepal.Width","Petal.Length","Petal.Width")

# 去掉标签列,只保留特征列
iris_features <- iris[1:4] # 另一种写法iris[-5]
# 构建一个函数,样本和邻居的点距
dist_eucl <- function(x1, x2) sqrt(sum((x1 - x2) ^ 2)) 
# 计算new_sample到iris_features150个样本的距离
distances <- apply(iris_features,1,function(x) dist_eucl(x,new_sample))
# 对计算后的距离排序
# index.return = T把计算的距离和数据框的行号对应起来
distances_sorted <- sort(distances,index.return = T) 
str(distances_sorted)

# 选取距离最近的前5个点
nn_5 <- iris[distances_sorted$ix[1:5],] 
# ix包含了对应观测数据的行号


5个最邻近的新样本中有4个属于versicolor品种,1个是virginica品种。通常采用投票作为筛选类别的预测。因此,认为测试样本点属于versicolor。
注:通常k取奇数个,方便投票。


1.2 模型的类别

有监督学习( supervised learning )
◾ Iris数据集,4个特征和1个输出量
◾ 对输入点的预测是通过组合该点的一小批邻近点的输出变量值来进行的
无监督学习(unsupervised learning)
◾ 若iris数据集的无监督版本只包含4个特征,没有品种的输出变量
◾ 查看数据并根据他们之间的相似度对观测数据进行分组
◾ 进行聚类(clustering)然后进行预测
半监督学习(semi-superivsed)
◾ 利用一部分包含输出量值而其他部分完全无标签的数据来粗略地训练模型,然后引入无标签数据,用训练的模型对其产生预测的标签
强化学习(reinforcement learning)
◾ 无法获得输出变量,但提供了其他和输出变量直接相关的信息根据棋谱,对赢得象棋比赛的下一步的走法进行预测

监督学习vs无监督学习
监督学习:输入{𝒙𝒊, 𝒚𝒊},学习𝑦 = 𝑓(𝑥; 𝜃)
分类: 𝑦是可分类的,例如:主题
回归: 𝑦是连续的,例如:温度,收入
排序: 𝑦是顺序的
结构预测:y是关于x的临时结构
无监督学习:输入{𝒙𝒊},学习𝑦 = 𝑓(𝑥; 𝜃)
概率密度估计: 𝑦是概率密度,𝑦 = 𝑝(𝑥; 𝜃)
异常监测: 𝑦是异常的

2. 预测建模的过程

定义模型的目标
收集数据
选取模型
数据预处理
特征工程和降维
训练和评估模型
重复尝试不同模型及模型的最终选择
部署模型

2.1数据预处理

◾探索性数据预处理
R函数中的,summary
图形可视化,ggplot2
psych包中,describe()函数
◾特征工程和降维
数值特征是完全互不相同的进位制来衡量的
输入特征进行变换(比例放缩)

2.1.1归一化和标准化


Caret包(Classification And REgression Training)
利用caret包的两步:
◾preProcess()函数:预处理的变换,训练数据集中
predict()函数:

  • 函数形式
    predict(object, newdata = NULL,
    type = c("link", "response", "terms"),
    se.fit = FALSE, dispersion = NULL, terms = NULL,na.action = na.pass, ...)
  • 参数介绍
    object: 是一个从glm继承的模型对象
    newdata: 是一个数据框,如果缺失就用训练数据进行预测,也即拟合值
    type: 表示预测种类。默认是归一化的线性预测;responses是归一化的响应变量。因此对于一个二分类模型,默认是log-odds (logit归一化的概率),然而type="response"给出的是预测概率。“terms”返回一个矩阵提供在线性预测下模型公式中每一项的拟合值。
    se.fit: 是一个bool值,表示标准误差是否需要。
    dispersion: GLM中用于计算标准化误差的离差。如果缺失,将会返回模型对象中的summary.
    terms: with type = "terms" by default all terms are returned. A character vector specifies which terms are to be returned.
    na.action: 表示如何对待newdata中的缺失数据,默认是将缺失值预测为NA.
library("caret")
# 去掉标签列,只保留特征列
iris_numeric <- iris[1:4]

# 区间归一化
pp_unit <- preProcess(iris_numeric, method =c("range")) 
iris_numeric_unit <- predict(pp_unit, iris_numeric)

# Z评分归一化
pp_zscore <- preProcess(iris_numeric,method=c("center","scale"))
iris_numeric_zscore <- predict(pp_zscore,iris_numeric)

# boxcox变换
pp_boxcox <- preProcess(iris_numeric,method=c("BoxCox"))
iris_numeric_boxcox <- predict(pp_boxcox,iris_numeric)
# 画图
# 归一化之前
ggplot(iris_numeric, aes(x=Sepal.Length)) + 
  geom_density(color="black",fill="black",alpha=0.4)+ 
  ggtitle("Unnormalized")

# 区间归一化
ggplot(iris_numeric_unit, aes(x=Sepal.Length)) + 
  geom_density(color="black",fill="black",alpha= 0.4)+ 
  ggtitle("Unit Interval Transformation")

# Z评分归一化
ggplot(iris_numeric_zscore, aes(x=Sepal.Length)) + 
  geom_density(color="black",fill="black",alpha= 0.4)+ 
  ggtitle("Z-Score Transformation")

# bocox变换
ggplot(iris_numeric_boxcox, aes(x=Sepal.Length)) + 
  geom_density(color="black",fill="black",alpha= 0.4)+ 
  ggtitle("Box Cox Transformation")
归一化之前

区间归一化

Z评分归一化

bocox变换
2.1.2寻找强相关性

◾caret包提供了findCorrelation()函数,可以把一个相关系数矩阵作为输入,还可以设定相关系数的绝对值指定一个阈值的cutoff参数。
◾它会返回一个(有可能长度为0)的向量。该向量会显示由于相关性需要从数据框中删除的列,cutoff的默认值是0.9。

# 查看各特征属性间的相关系数,并将相关系数赋iris_cor变量
cor(iris_numeric)

可见,Petal.Width和Petal.Length的相关系数高达0.96。

iris_cor <- cor(iris_numeric)
findCorrelation(iris_cor) 
findCorrelation(iris_cor,cutoff=0.99) 
findCorrelation(iris_cor,cutoff=0.80) 

由于findCorrelation默认是找到相关系数大于90%的特征,可以看到,返回结果为3,代表Petal.Length 与其他变量存在较强相关性,可以被删除。
如果设置cutoff=0.99,结论是0,找不到需要剔除的变量;如果设置cutoff=0.8,将发现要剔除的变量有两个。

2.1.3寻找共线性

两个变量强相关,并不能说明他们一定具有共线性。如要发现变量间的共线性,可以使用findLinearCombos函数。

#构造一个新的数据框,人为增加两列。
#Cmb是Sepal.Length和Petal.Width的线性表示,Cmb.N是在Cmb上再加一个随机数
new_iris <- iris_numeric
new_iris$Cmb <- 6.7*new_iris$Sepal.Length - 0.9*new_iris$Petal.Width
set.seed(68)
new_iris$Cmb.N <- new_iris$Cmb + rnorm(nrow(new_iris),sd=0.1)
options(digits = 4)  #为了显示方便,只显示四位数
head(new_iris, n = 3) 

Cmb是Sepal.Length和Petal.Width特征的完全线性组合,Cmb.N和Cmb相同,但是加上了某些均值为0,且标准差很小(0.1)的高斯噪声的特征。
利用findLinearCombos()函数,可以检测出严格的特征线性组合,不过有噪声就检测不到了。

# 使用findLinearCombos函数来查找共线性
findLinearCombos(new_iris)


可以看到,建议删除第5个特征(Cmb),因为它是1和4的线性组合。

2.1.4寻找零方差变量

零方差变量就是值为一个常数的变量,所有观测数据在该特征属性上都是一个常数,或者绝大部份都是一个常数,导致数据在抽样过程中没有变化。这样的数据会影响模型的准确性。
在R中,使用nearZeroVar函数来检测零方差变量或则近似零方差变量。

# 构造一个新的数据框,人为增加两列。
# ZV列是一个常数6.5(零方差变量)
# Yellow除了第一行是True,其他都是False,所以是一个近似零方差变量
newer_iris <- iris_numeric
newer_iris$ZV <- 6.5
newer_iris$Yellow <- ifelse(rownames(newer_iris)==1,T,F)

nearZeroVar(newer_iris)
nearZeroVar(newer_iris, saveMetrics=T)

freqRatio频数比率,percentUnique 唯一值比率,zeroVar表示零方差变量,nzv表示近似零方差变量。

2.1.5特征降维——PCA

R语言有非常多的途径做主成分分析,比如自带的princomp()和prcomp()函数,还有psych包的principal()函数、gmodels包的fast.prcomp函数、FactoMineR包的PCA()函数、ade4包的dudi.pca()函数、amap包的acp()函数等。
主成分分析法(principal component analysis)
◾主成分分析经常用减少数据集的维数,同时保持数据集的对方差贡献最大的特征
◾原理:
设法将原来变量重新组合成一组新的相互无关的几个综合变量,同时根据实际需要从中可以取出几个较少的综合变量尽可能多地反映原来变量的信息的统计方法,也是处理降维的一种方法。
◾选取前两个主成分或其中某两个主成分,由图形可直观地看出各个输入在主分量中的地位。

# thresh表示PCA的累积方差比达到95%
pp_pca <- preProcess(iris_numeric, method=c("BoxCox", "center", "scale", "pca"), thresh=0.95)
iris_numeric_pca <- predict(pp_pca, iris_numeric)
#保留两位数字
options(digits=2)  
# 使用pp_pca$rotation来查看变量的有效载荷(权值),或特征向量、回归系数
pp_pca$rotation

Iris花的数值特征的前两个主成分就包含了数据集里面95%的变异。



主成分权值,也就是说:
PC1=𝟎. 𝟓𝟐 ∗ Se𝒑𝒂𝒍. 𝑳𝒆𝒏𝒈𝒕𝒉 − 𝟎. 𝟐𝟕 ∗ 𝑺𝒆𝒑𝒂𝒍.𝑾𝒊𝒅𝒕𝒉 + 𝟎. 𝟓𝟖 ∗ 𝑷𝒆𝒕𝒂𝒍. 𝑳𝒆𝒏𝒈𝒕𝒉 + 𝟎. 𝟓𝟕 ∗ 𝑷𝒆𝒕𝒂𝒍. 𝑾𝒊𝒅𝒕

pp_pca_full <- preProcess(iris_numeric, method=c( "BoxCox", "center", "scale","pca"),pcaComp=4) 
# pcaComp=4,用4个PCA维度
iris_pca_full <- predict(pp_pca_full, iris_numeric)
 #150行数据形成PCA下的新坐标
head(iris_pca_full,3)

Iris花的数值特征的前两个主成分就包含了数据集里面95%的变异:
第一个PCA是一条拟合原始数据最近似的直线
第二个PCA是拟合第一条直线拟合原始数据的误差
第三个PCA是拟合第一个和第二个PCA拟合后剩下的误差

#计算PC1、PC2、PC3、PC4 四个维度下的方差
pp_pca_var <- apply(iris_pca_full,2,var) 
pp_pca_var 
#计算每一列的百分占比;计算每一列的累计百分占比
iris_pca_var <- data.frame(Variance = round(100*pp_pca_var/sum(pp_pca_var),2), CumulativeVariance = round(100*cumsum(pp_pca_var) / sum(pp_pca_var),2)) 
iris_pca_var

计算得出来PC1新坐标的占比是73%,PC1和PC2占比96%。也就是PC1和pC2两个新的维度可以解释原数据96%。

# 画主成分的陡坡图
p <- ggplot(data=iris_pca_var) 
p <- p + aes(x=rownames(iris_pca_var), y=Variance, group=1) #x轴是PC1-PC4
p <- p + geom_line() #连线
p <- p + geom_point() #点描出来
p <- p + xlab("Principal Component") #横坐标是PC1-PC4
p <- p + ylab ("Percentage of Original Variance Captured") #纵坐标是PC1-PC4对应方差
p <- p +  ggtitle("Scree Plot for Iris PCA") 
p <- p + geom_text(aes(label = paste(Variance, "% (", CumulativeVariance, "%)", sep = ""), parse = F),hjust=0, vjust=-0.5)
主成分的陡坡图

2.2训练模型

训练集(training set)和测试集(test set)
◾通常测试集占15%~30%
◾如果数据较少划分为85%/15%
◾如果数据较多划分为70%/30%
过拟合(overfitting)
◾在训练集中表现良好,但是在测试集表现较差
◾对观测数据拟合得过于严密,无法对未知数据广义化
◾换言之,该模型在训练数据集中提取了不真实的细节和变异,而这些细节和变异不具有代表性
偏误-方差权衡(bias-variance tradeoff)
◾训练和测试数据性能不一致的来源还包括模型偏误和方差
◾统计模型的方差:如果使用另一个训练集,模型预测出来的函数会有多大的变化。我们希望方差越小越好。

其中,caret包中用 createDataPartition()函数 创建抽样向量。

在非监督学习中,使用的是sample函数来进行随机抽样。
sample(x, size, replace = FALSE, prob = NULL)
x表示要抽样的向量,size表示要抽样的次数,replace代表是否是放回抽样,默认是不放回抽样。prob也是一个向量,用于标示待抽样向量中每个元素的权重。

# 随机抽样,该函数将会随机抽取80%的数据。list=false表示返回是一个矩阵而不是列表。
set.seed(2412)
# 按照种类进行索引,提取80%作为训练集
iris_sampling_vector <- createDataPartition(iris$Species, p = 0.8, list = FALSE) 
iris_sampling_vector

# 获取训练集数据
iris_train <- iris_numeric[iris_sampling_vector,]
iris_train_z <- iris_numeric_zscore[iris_sampling_vector,] # Z评分归一化,选取80%训练集
iris_train_pca <- iris_numeric_pca[iris_sampling_vector,] # PCA降维后,选取80%作为训练集
iris_train_labels <- iris$Species[iris_sampling_vector]

# 获取测试集数据,去掉训练集所在的行
iris_test <- iris_numeric[-iris_sampling_vector,] #iris_test_z <- iris_numeric_zscore[-iris_sampling_vector,] # Z评分归一化后测试集
iris_test_pca <- iris_numeric_pca[-iris_sampling_vector,] # PCA降维后测试集
iris_test_labels <- iris$Species[-iris_sampling_vector]

# 找5个邻居模式下,用原始数据集iris_train训练得到knn_model
knn_model <- knn3(iris_train, iris_train_labels, k = 5) 
# 找5个邻居模式下,用Z评分归一化数据集iris_train_z训练得到knn_model
knn_model_z <- knn3(iris_train_z, iris_train_labels, k = 5)
# 找5个邻居模式下,用PCA降维数据集iris_train_pca训练得到knn_model
knn_model_pca <- knn3(iris_train_pca, iris_train_labels, k = 5)
# 画图
library(class)
require(class)

layout(matrix(1:4,2,2, byrow=T)) #4个图,2行2列显示

for (k in c(1,5,10,15)) {
    x <- iris_train_pca
    px1 <- seq(from=-3.5, to=3.5, by=0.1)
    px2 <- seq(from=-3.5, to=3.5, by=0.1)
    g <- iris_train_labels
    xnew <- as.matrix(expand.grid(px1, px2))
    mod15 <- knn(x, xnew, g, k=k, prob=F)
    prob <- ifelse(mod15=="setosa", 1.0, ifelse(mod15=="virginica", 0.6, 0.1))
    prob15 <- matrix(prob, length(px1), length(px2))
    par(mar=rep(2,4))
    contour(px1, px2, prob15, levels=c(0.2), labels="", xlab="", ylab="", main=paste(k,"NN on Iris PCA",sep=""), axes=FALSE)
    train_pch <- ifelse(g=="setosa", 15, ifelse(g=="virginica", 16, 17))
    grid_pch <- ifelse(mod15=="setosa", 15, ifelse(mod15=="virginica", 16 , 17))
    points(x, pch=train_pch, cex=1.2, col="gray40")
    gd <- expand.grid(x=px1, y=px2)
    points(gd, pch=".", cex=0.4)
    box()
}

2.3评估模型



卡巴统计量是被设计用来抵消随机因素的,它的取值空间为[-1, 1],其中1代表完全精确,-1代表完全不精确,而0是在精确度恰好等于随机猜测法的结果时出现的。
它主要用于一致性检验,用来度量两个被观测对象的一致程度。
在R中,使用caret包中的postResample函数来计算准确率和卡巴系数。

knn_predictions <- predict(knn_model, iris_test, type = "class")
postResample(knn_predictions, iris_test_labels)
# 还可以使用table来显示混淆矩阵(交叉列联表)
table(knn_predictions,iris_test_labels)


评估二元分类问题
第一类错误(Type I error):假阳性
第二类错误(Type II error):假阴性


在R中,其实可以使用caret包的knn3函数来方便的进行分类预测。下面是完整代码:

#剔除标签项,保留特征列
iris_numeric <- iris[1:4]
#数据预处理,区间归一化
pp_unit <- preProcess(iris_numeric, method = c("range"))
iris_numeric_unit <- predict(pp_unit, iris_numeric)

#随机抽样。划分训练集和测试集
iris_sampling_vector <- createDataPartition(iris$Species, p = 0.8, list = FALSE)
#获取训练集数据
iris_train <- iris_numeric_unit[iris_sampling_vector,]
#获取训练集分类数据
iris_train_labels <- iris$Species[iris_sampling_vector]
#获取测试集数据
iris_test <- iris_numeric_unit[-iris_sampling_vector,]
#获得测试集分类数据
iris_test_labels <- iris$Species[-iris_sampling_vector]

#训练模型
knn_model <- knn3(iris_train, iris_train_labels, k=5)

#根据模型对测试集进行预测,评分类型为分类或者投票
knn_predictions <- predict(knn_model, iris_test, type = "class")

#计算准确率
postResample(knn_predictions, iris_test_labels)

#显示混淆矩阵
table(knn_predictions,iris_test_labels)
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 204,793评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,567评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,342评论 0 338
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,825评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,814评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,680评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,033评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,687评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 42,175评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,668评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,775评论 1 332
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,419评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,020评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,978评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,206评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,092评论 2 351
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,510评论 2 343

推荐阅读更多精彩内容