搜索此博客

2017年6月27日星期二

Rmd语法

文本输出
代码高亮(highlight=TRUE),增强可读性,有无数的高亮主题可选,仅适用于LaTeX和HTML输出,MD文档在转为HTML文档之后可以用专门的JavaScript库去高亮代码
代码重排(tidy=TRUE),对那些不注意代码格式的人来说很有用,再乱的代码,到了这里也会变得相对整齐,本功能由formatR包支持
执行或不执行代码(eval=TRUE/FALSE),不执行的代码段将被跳过,原样输出源代码
显示/隐藏源代码(echo=TRUE/FALSE),甚至精确控制显示哪几段代码(echo取数值)
显示/隐藏普通文本输出或将文本输出以原样形式输出(results='markup', 'hide', 'asis')
显示/隐藏警告文本(warning=TRUE/FALSE)、错误消息(error)和普通消息(message)
显示/隐藏整个代码段的输出(include=TRUE/FALSE),比如我们可能想运行代码,但不把结果写入输出中
child选项 加选项child='文件名.Rmd'可以调入另一个.Rmd文件的内容。 如果有多个.Rmd文件依赖于相同的代码,可以用这样的方法。
collapse选项 一个代码块的代码、输出通常被分解为多个原样文本块中, 如果一个代码块希望所有的代码、输出都写到同一个原样文本块中, 加选项collapse=TRUE。 
results选项 用选项results=选择文本型结果的类型。 取值有:
markup, 这是缺省选项, 会把文本型结果变成HTML的原样文本格式。
hide, 运行了代码后不显示运行结果。
hold, 一个代码块所有的代码都显示完, 才显示所有的结果。
asis, 文本型输出直接进入到HTML文件中, 这需要R代码直接生成HTML标签, knitr包的kable()函数可以把数据框转换为HTML代码的表格。

输出表格
knitr包提供了一个 kable() 函数可以用来把数据框或矩阵转化成有格式的表格, 支持HTML、docx、LaTeX等格式。 因为kable()的结果已经是HTML格式的, 为了将其结果直接插入到HTML内容中, 需要使用代码段选项results='asis'。如:
x <- 1:10; y <- x^2; lmr <- lm(y ~ x)
co <- summary(lmr)$coefficients
print(co)
print(knitr::kable(co))
扩展包xtable提供了一个xtable()函数, 也可以用来生成HTML格式和LaTeX格式的表格, 但是需要指定要输出的格式。 xtable对比较多的R数据类型和输出类型提供了表格式显示功能, 包括矩阵、数据框、回归分析结果、方差分析结果、主成分分析结果、 若干分析结果的summary结果等。 例如,上面的回归结果用xtable()函数显示如:
print(xtable::xtable(lmr), type='html')
意这里指定了输出为HTML类型。 如果将本文件转化为docx, xtable的结果不可用。

图形控制
fig.path用来设置图形输出的路径,对大型报告来说,我们可能希望把各种乱七八糟的文件归类管理,所以R图形文件可以用这个选项写入单独的文件夹
dev设置用哪种图形设备记录图形,自带二十多种常见的图形设备,如PDF、PNG甚至tikz,具体取值参见网站上的文档
fig.width和fig.width设置图形文件本身的宽高尺寸(单位英寸)
out.width和out.height设置图片在输出文档中的宽高,这是相对Sweave来说的新选项,也是我等了很久没等到,最终刺激我自己写包的原因之一(这两个选项太容易实现了),很多用户看到这两个选项都很欣喜,一个小问题,困扰多少英雄好汉
fig.keep设置保留图形的方式,我们可以完全不保留,也可以只保留高层作图函数生成的图形,也可以保留低层作图函数产生的图形
fig.show设置图形显示的方式,可以跟在作图代码后面即刻显示出来,也可以等到代码段运行完毕之后再把该段中所有图形一气儿显示出来,也可以把所有图形显示为动画

参考文献
http://www.math.pku.edu.cn/teachers/lidf/docs/Rbook/rmarkdown.html
https://github.com/yihui/r-ninja/blob/master/11-auto-report.Rmd

2017年6月21日星期三

kerasR包学习

安装

devtools::install_github("statsmaths/kerasR")

例子(波士顿住房数据)

Keras建立一个模型,首先构建一个空的Sequential模型。
library(kerasR)
## successfully loaded keras
mod <- Sequential()
Sequential的结果与kerasR提供的大多数功能一样,是一个python.builtin.object。这种从reticulate包定义的对象类型可以直接访问底层python类暴露的所有方法和属性。要访问这些,需要使用$运算符后跟方法名称,通过调用add方法添加层。此函数作为另一个python.builtin.object的输入,通常构造为另一个kerasR函数的输出。例如,要在我们的模型中添加一个dense 层,我们执行以下操作:
mod$add(Dense(units = 50, input_shape = 13))
我们现在已经添加了一个具有200个神经元的致密层。 第一层必须包含input_shape的规范,给出输入数据的维度。 这里我们设置输入变量的数量等于13.接下来在模型中,我们将添加relu激活到模型中:
mod$add(Activation("relu"))
现在,我们添加一个只有一个神经元的密集层作为输出层:
mod$add(Dense(units = 1))
一旦模型被完全定义,我们必须在拟合它的参数或使用它进行预测之前进行编译。编译模型可以使用方法compile来完成,但是它的一些可选参数在从R类型转换时会引起麻烦,所以我们提供了一个自定义的包装器keras_compile。 至少我们需要指定损失函数和优化器。损失可以用一个字符串来指定,但是我们将把另一个kerasR函数的输出作为优化器传递。 这里我们使用RMSprop优化器,因为它通常具有相当好的性能:
keras_compile(mod, loss = 'mse', optimizer = RMSprop())
现在我们可以从一些训练数据中将模型中的权重进行拟合,但是我们还没有从中训练的任何数据。我们使用load_boston_housing函数加载数据。 我们提供几种数据加载功能作为软件包的一部分,并且以相同的格式返回数据。这个例子中,矩阵进行归一化是有意义。虽然我们可以使用R函数进行归一化,在这里我们使用keras特定的函数进行归一化。 归一化的一个好处是不同纬度的测量可以统一,这是卷积和循环神经网络中的有用特征。
boston <- load_boston_housing()
X_train <- normalize(boston$X_train)
Y_train <- boston$Y_train
X_test <- normalize(boston$X_test)
Y_test <- boston$Y_test
现在,我们使用keras_fit从数据中拟合模型。与编译一样,有一种直接的方法,但是您可能会遇到直接调用它的数据类型转换问题。 相反,直接使用包装函数功能比较容易(如果您自己运行,您将看到Keras为跟踪模型的拟合提供了非常好的详细输出):
keras_fit(mod, X_train, Y_train,batch_size = 32, epochs = 200,verbose = 1, validation_split = 0.1)
请注意,该模型在这里并不特别好,这可能是由于在小数据集上产生了过度拟合。
pred <- keras_predict(mod, X_test)
sd(as.numeric(pred) - Y_test) / sd(Y_test)
## [1] 0.8261666

另外还有几个例子包含在包装Vignette RKeras Deep Learning Library中。

R interface to Keras

R interface to Keras
有多个版本的python时需要选择一个
sudo pip3 install virtualenv
#首先罗列出所有可用的 python 替代版本信息
update-alternatives --list python
#update-alternatives: error: no alternatives for python
如果出现以上所示的错误信息,则表示 Python 的替代版本尚未被 update-alternatives 命令识别

#--install 选项使用了多个参数用于创建符号链接。最后一个参数指定了此选项的优先级,如果我们没有手动来设置替代选项,那么具有最高优先级的选项就会被选中。
sudo update-alternatives --install /usr/bin/python python /usr/bin/python2.7 1
sudo update-alternatives --install /usr/bin/python python /usr/bin/python3.5 2
# update-alternatives --config python
python --version

安装

devtools::install_github("rstudio/keras")
library(keras)
install_tensorflow()
#安装GPU版本的tensorflow install_tensorflow(gpu=TRUE),https://rstudio.github.io/keras/

快速入门

Keras的核心数据结构是一个模型,是组织层次的一种方式。 最简单的模型是Sequential模型,一个线性层叠层。 对于更复杂的体系结构,您应该使用Keras功能API,这允许构建层的任意图形。 定义一个顺序模型:
library(keras)
model <- keras_model_sequential()
model %>%
layer_dense(units = 512, activation = 'relu', input_shape = c(784)) %>%
layer_dropout(rate = 0.2) %>%
layer_dense(units = 512, activation = 'relu') %>%
layer_dropout(rate = 0.2) %>%
layer_dense(units = 10, activation = 'softmax')
#使用summary()函数打印模型结构的摘要:
summary(model)
## Model
## ___________________________________________________________________________
## Layer (type) Output Shape Param #
## ===========================================================================
## dense_1 (Dense) (None, 512) 401920
## ___________________________________________________________________________
## dropout_1 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_2 (Dense) (None, 512) 262656
## ___________________________________________________________________________
## dropout_2 (Dropout) (None, 512) 0
## ___________________________________________________________________________
## dense_3 (Dense) (None, 10) 5130
## ===========================================================================
## Total params: 669,706
## Trainable params: 669,706
## Non-trainable params: 0
## ___________________________________________________________________________
##
##
采用适当的损失函数,优化器和矩阵创建模型:
model %>% compile(
loss = 'categorical_crossentropy',
optimizer = optimizer_sgd(lr = 0.02),
metrics = c('accuracy')
)
现在可以批量地迭代训练数据(x_trainy_train是矩阵):
history <- model %>% fit(
x_train, y_train,
epochs = 20, batch_size = 128,
validation_split = 0.2
)
绘制来自训练数据的损失和准确性指标:
plot(history)
评估测试数据的性能:
loss_and_metrics <- model %>% evaluate(x_test, y_test, batch_size = 128)
生成关于新数据的预测:
classes <- model %>% predict(x_test, batch_size = 128)

R中使用MLPMNIST手写数字进行分类

#https://www.analyticsvidhya.com/blog/2017/06/getting-started-with-deep-learning-using-keras-in-r/
library(keras)
#加载mnist数据集
data<-dataset_mnist()

#拆分训练集和测试集
train_x<-data$train$x
train_y<-data$train$y
test_x<-data$test$x
test_y<-data$test$y

rm(data)
2D数组转换为1D数组,送到MLP中,并进行归一化
train_x <- array(train_x, dim = c(dim(train_x)[1], prod(dim(train_x)[-1]))) / 255
test_x <- array(test_x, dim = c(dim(test_x)[1], prod(dim(test_x)[-1]))) / 255
使用keras内置函数将目标变量转换为one-hot编码
train_y<-to_categorical(train_y,10)
test_y<-to_categorical(test_y,10)
定义keras sequential model
model <- keras_model_sequential()
使用1个输入层[784个神经元]1个隐层[784个神经元],具有dropout 0.41个输出层[10个神经元]定义模型
model %>%
layer_dense(units = 784, input_shape = 784) %>%
layer_dropout(rate=0.4)%>%
layer_activation(activation = 'relu') %>%
layer_dense(units = 10) %>%
layer_activation(activation = 'softmax')
metric = precisionoptimiser作为adam编译定义的模型
model %>% compile(
loss = 'categorical_crossentropy',
optimizer = 'adam',
metrics = c('accuracy')
)
模型在训练数据集上拟合
model %>% fit(train_x, train_y, epochs = 100, batch_size = 128)
交叉验证数据集评估模型
loss_and_metrics <- model %>% evaluate(test_x, test_y, batch_size = 128)

上述代码的训练集准确率为99.14,验证准确率为96.89。 代码在i5处理器上运行,在单个epoch占用13.5秒,而在TITANx GPU上,验证精度为98.44,平均epoch2秒。

2017年6月20日星期二

formatR包学习

formatR很多种用法,这里只展示三种。
1.载入包后,直接在R的命令行窗口输入tidy_app() ,该命令就会调用一个shiny程序,自动在你的浏览器里打开一个界面模式下的格式化工具。
library(formatR)
#tidy_app()
2.Addins--Reformat R Code 该命令会对当前.R的文件自动排版,不支持.Rmd文件。

3.先复制你需要格式化的代码,在R命令行窗口输入tidy_source(),该函数就会自动输出格式化后的代码。

timekit包学习

http://rpubs.com/xuefliang/286234

2017年6月13日星期二

Ocserv 搭建 Cisco AnyConnect VPN服务端 一键脚本

  1. wget -N --no-check-certificate https://softs.pw/Bash/ocserv.sh && chmod +x ocserv.sh && bash ocserv.sh

  2. # 如果上面这个脚本无法下载,尝试使用备用下载:
  3. wget -N --no-check-certificate https://raw.githubusercontent.com/ToyoDAdoubi/doubi/master/ocserv.sh && chmod +x ocserv.sh && bash ocserv.sh

2017年6月11日星期日

lubridate包

R语言的基础包中提供了两种类型的时间数据,一类是Date日期数据,它不包括时间和时区信息,另一类是POSIXct/POSIXlt类型数据,其中包括了日期、时间和时区信息。一般来讲,R语言中建立时序数据是通过字符型转化而来,但由于时序数据形式多样,而且R中存贮格式也是五花八门,例如Date/ts/xts/zoo/tis/fts等等。用户很容易被一系列的数据格式所迷惑,所以时序数据的转化和操作并不是非常方便。所幸的是,我们有了lubridate包。lubridate包主要有两类函数,一类是处理时点数据(time instants),另一类是处理时段数据(time spans
时点类函数,它包括了解析、抽取、修改。
library(lubridate)
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
##
## date
# 从字符型数据解析时间,会自动识别各种分隔符
x <- ymd('2010-04-08')
# 观察x日期是一年中的第几天
yday(x)
## [1] 98
# 修改x日期中的月份为5
month(x) <- 5
时段类函数,它可以处理三类对象,分别是: interval:最简单的时段对象,它由两个时点数据构成。 duration:去除了时间两端的信息,纯粹以秒为单位计算时段的长度,不考虑闰年和闰秒,它同时也兼容基本包中的difftime类型对象。 period:以较长的时钟周期来计算时段长度,它考虑了闰年和闰秒,适用于长期的时间计算。以2012年为例,duration计算的一年是标准不变的365天,而period计算的一年就会变成366天。 有了时点和时段数据,就可以进行各种计算了。
# 从两个时点生成一个interval时段数据
y <- interval(x,now())
# interval格式转为duration格式
z <- as.duration(as.integer(as.difftime(y,units = "secs")))
# 时点+时段生成一个新的时点
now() + as.duration(z)
## [1] "2024-07-16 11:58:30 CST"
# 10天后的时间数据
now() + ddays(10)
## [1] "2017-06-21 21:59:15 CST"
# 下面的两条语句很容易看出durationperiod的区别,dyears(1)表示duration对象的一年,它永远是365天。而years(1)表示period对象的一年,它识别出2012是闰年,它有366天,所以得到正确的时点。
ymd('20120101') + dyears(1)
## [1] "2012-12-31"
ymd('20120101') + years(1)
## [1] "2013-01-01"
为了处理时区信息,lubridate包提供了三个函数: tz:提取时间数据的时区 with_tz:将时间数据转换为另一个时区的同一时间 force_tz:将时间数据的时区强制转换为另一个时区
# 输入欧洲杯决赛在乌克兰的开场时间,再转为北京时间
eurotime <- ymd_hms('2012-07-01 21:45:00',tz='EET')
with_tz(eurotime,tzone='asia/shanghai')
## [1] "2012-07-01 18:45:00 asia"
最后来玩玩股票指数作为结束吧,在金融市场中谣传着一种日历效应。即在一周的第一天或者一年的第一个月份,股票会出现不错的上涨。让我们用热图来观察一下。首先是获取上证指数数据,然后根据不同的月份和星期数,将收益率汇集到不同的组中。将该组收益率的中位数映射到图形颜色上去。可以从下图看到,似乎并不存在明显周一效应或是一月效应。
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: TTR
## Version 0.4-0 included new data defaults. See ?getSymbols.
library(ggplot2)
library(lubridate)
# 读取上证指数历史数据
getSymbols('^SSEC',src='yahoo',from = '1997-01-01')
## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.
##
## WARNING: There have been significant changes to Yahoo Finance data.
## Please see the Warning section of '?getSymbols.yahoo' for details.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.yahoo.warning"=FALSE).
## Warning: ^SSEC contains missing values. Some functions will not work if
## objects contain missing values in the middle of the series. Consider using
## na.omit(), na.approx(), na.fill(), etc to remove or replace them.
## [1] "SSEC"
time <- ymd(as.character(index(SSEC)))
open <- as.numeric(Op(SSEC))
high <- as.numeric(Hi(SSEC))
low <- as.numeric(Lo(SSEC))
close <- as.numeric(Cl(SSEC))
volume <- as.numeric(Vo(SSEC))
# 根据收盘和开盘计算当日收益率
profit <- (close-open)/open
# 提取时间数据中的周数和月份
wday <- wday(time)-1
mday <-month(time)
data <- data.frame(time,wday,mday,profit)
p <- ggplot(data,aes(factor(mday),factor(wday),z=profit))
# 收益率的热图,图中颜色越浅,表示汇集到这个组中的收益率中位数越高
p +stat_summary_2d(fun=function(x) median(x))+ labs(x='月份',y='星期')
## Warning: Removed 172 rows containing non-finite values (stat_summary2d).