灰色模型的R代码

最近帮朋友写了一个灰色模型GM(1,1)的R实现，参考网上现有的matlab代码，比较容易就可以弄出来。下面是具体过程，主函数是GM()，建立的模型是一个S3类，搭配两个自定义的泛型函数print和plot可以得到结果输出和图形。

# 本代码用来计算GM(1,1)灰度模型
# 输入：x 原始数据，adv为外推预测步长
# 输出：actual 原始数据，fit 拟合数据,degree 拟合精度，
#       C 后验差比值，P 小概率误差，predict 外推预测值
 
 
# 主函数
GM <- function(x,adv=0) {
    x0 <- x
    k = length(x0)
    # AGO
    x1 = cumsum(x0)
    # construct matrix B & Y
    B = cbind(-0.5*(x1[-k]+x1[-1]),rep(1,times=k-1))
    Y = x0[-1]
    # compute BtB...
    BtB = t(B)%*%B
    BtB.inv = solve(BtB)
    BtY = t(B)%*%Y
 
    # estimate
    alpha = BtB.inv%*%BtY
 
    # 建立预测模型
 
    predict <- function(k) {
        y = (x0[1] - alpha[2]/alpha[1])*exp(-alpha[1]*k)+alpha[2]/alpha[1]
        return(y)
    }
    pre <- sapply(X=0:(k-1),FUN=predict)
    prediff <- c(pre[1],diff(pre))
    # 计算残差
    error <- round(abs(prediff-x0),digits=6)
    emax <- max(error)
    emin <- min(error)
    # 模型评价
    incidence <- function(x) {
         return((emin + 0.5*emax)/(x+0.5*emax))
    }
    degree <- mean(sapply(error,incidence))
 
    s1 <- sqrt(sum((x0-mean(x0))^2)/5)
    s2 <- sqrt(sum((error-mean(error))^2)/5)
 
    C <- s2/s1
 
    e <- abs(error-mean(error))
    p <- length(e<0.6745)/length(e)
 
    result <- list(actual = x0,
                   fit = prediff,
                   degree = degree,
                   C = C,
                   P = p)
 
    # 外推预测第k+adv项
    if (adv > 0) {
        pre.adv <- predict(k+adv-1)-predict(k+adv-2)
 
        result$predict <- pre.adv
     }
    class(result) <- 'GM1.1'
    return(result)
}
 
# 增加对应gm1.1类的泛型函数 print & plot
print.GM1.1 <- function(mod){
    cat('the result of GM(1,1)\n')
    cat('Actual Data:', '\n',mod$actual ,'\n')
    cat('Fit Data:', '\n',round(mod$fit,2) ,'\n')
    cat('Degree:', round(mod$degree,3) ,'\n')
    cat('C:', round(mod$C,3) ,'\n')
    cat('P:', round(mod$P,3) ,'\n')
    if (!is.null(mod$predict)) {
        cat('Predict Data:', round(mod$predict,2), '\n')
    }
}
 
plot.GM1.1 <- function(mod,adv=5){
    prex <- numeric(adv)
    x <- mod$actual
    for (k in 1:adv){
        prex[k] <- GM(x,k)$predict    
    }
 
    value = c(x,prex)
 
    res <- data.frame(index = 1:length(value),
                      value = value,
                      type = factor(c(rep(1,length(x)),rep(2,length(prex)))))
    library(ggplot2)
    p <- ggplot(res,aes(x=index,y= value))
    p + geom_point(aes(color=type),size=3)+ 
        geom_path(linetype=2) +
        theme_bw()
}
 
 
# 原始数据
x = c(26.7,31.5,32.8,34.1,35.8,37.5)
 
# 预测第7项
res <- GM(x,1)
print(res)
plot(res,3)

那些奇葩的R函数

看别人的代码会遇到一些奇葩的函数，一般的教程上很少提到，但却有很好的用处，这类函数基本上分布在base以及utils包中，下面将它们略为归纳一下，以备后用。

1，文件执行：
在用R生成一个PDF文档后，如果想去打开它，你可能会在文件夹里找到再点开。再或者我们想调用系统中的其它程序来做点事情，可能要打开cmd敲点命令。实际上这都可以在R内部完成。举例来说用pandoc转换na.md成docx再打开它。
system('pandoc d:\\rspace\\na.md -o d:\\rspace\\na.docx')
shell.exec('d:\\rspace\\na.docx')

2，网络浏览：
browseURL：浏览某个指定的网页
download.file：下载网络文件到本地

3，文件操作
dir.create：新建一个文件夹
list.dirs：显示目录下的文件夹
list.files：显示目录下的文档
file.create：文档创建
file.exists：判断文档是否存在
file.remove：文档删除
file.rename：重命名
file.append：文档添加
file.copy：文档复制
file.symlink(from, to)
file.show：显示文档内容
file.info：显示文档信息
file.edit：编辑文档
zip： 压缩文件
unzip： 解压缩文件

4，运算进度条
在一个大循环运算时，如果可以看到目前的进度是比较方便的，txtProgressBar和setTxtProgressBar函数可以帮助做到这一点，下面是内置的一个小例子：

testit <- ...="..." function="function" p="p" x="sort(runif(20)),">{
    pb <- p="p" txtprogressbar="txtprogressbar">    for(i in c(0, x, 1)) {Sys.sleep(0.5); setTxtProgressBar(pb, i)}
    Sys.sleep(1)
    close(pb)
}
testit()

笨办法学R编程（6）

有时候用R来解一些Project Euler的题目会非常简单，今天就来三题连解（6、7、8）。题目就不再这里复述了，可以查看官方网站。用到函数和表达式大部分在前面都已经熟悉过了，不过还是会接触到一些新的函数。废话不多说直接上代码吧。

# Euler 6  
x <- 1:100
sum(x)^2 - sum(x^2)
 
# Euler 7 
 
n <- 0
i <- 1
m <- rep(0,10001)
while (n <10001) {
    if (findprime(i)) {
        n <- n +1 
        m[n] <- i}
    i <- i + 1
}
m[10001]
 
# 预备练习，熟悉一些字符串操作函数
text <- c('hello','world','I','love','code')
gsub('o',' ',text)
gsub('o','*',text)
gsub('o','',text)
 
(temp1 <-paste(text,collapse=' '))
paste(text,collapse='*')
paste(text,collapse='')
 
(temp2 <- strsplit(temp1," "))
class(temp2)
(temp3 <- unlist(strsplit(temp1," ")))
class(temp3)
 
# Euler 8 
 
web <- 'http://projecteuler.net/problem=8'
# 用readLines函数来抓取网页
raw <- readLines(web)
raw <- raw[54:73]
# 删除多余字符串
data <- gsub('<br />','',raw)
# 粘合成一个字符串
num <- paste(data,collapse='')
# 分割后转为数值向量
temp <- as.numeric(unlist(strsplit(num,'')))
 
p <- numeric()
for ( i in 1:(1000-4)) {
    p[i] <- prod(temp[i:(i+4)])
}
max(p)

第六题用R的向量化计算非常简单，第七题需要用到第二题中建立的findprime函数。第八题的解决是用字符串方法，先将那一长串数字作为字符串切开，再转为数值型向量，最后用循环求乘积。为了偷懒，是直接抓取的网页，没有输入那个长长的数字。

笨办法学R编程（5）

随着教程推进，基本的语法都接触得差不多了。当要解决某个具体问题时，只需要考虑用什么样的算法来整合运用这些函数和表达式。今天来解决Project Euler的第五个问题，该问题可以用很笨的暴力搜索法子来作，但是更聪明的作法是采用质因子分解的思路。即任何一个合数都可以分解为质数的乘积。为了完成这个题目，还需要学习一点点矩阵，以及和sapply函数相似的另一个函数apply。

# 预备练习
mat <- matrix(1:12,ncol=4)
print(mat)
t(mat)
colnames(mat) <- c('one','two','three','four')
rownames(mat) <- c('a','b','c')
print(mat)
apply(mat,1,sum)
apply(mat,2,sum)
sum(apply(mat,2,sum))
prod(apply(mat,2,sum))

笨办法学R编程（4）

看到各位对“笨办法系列”的东西还比较感兴趣，我也很乐意继续写下去。今天的示例将会用到数据框（data.frame）这种数据类型，并学习如何组合计算两个向量，以及如何排序。我们将用所学的东西来解决Project Euler的第四个问题，就是找出一个集合中最大的回文数。回文数是指一个像1534351这样“对称”的数，如果将这个数的数字按相反的顺序重新排列后，所得到的数和原来的数一样。开始啦！

# 预备练习
 
x <- y <- 1:9
data <- expand.grid(x=x,y=y)
print(data)
z <- data$x * data$y
# 一个九九乘法表
z <- matrix(z,ncol=9)
 
set.seed(1)
x <- round(runif(10),2)
print(x)
order(x)
x[order(x)[1]]
which.min(x)
x[which.min(x)]
x[order(x)]

笨办法学R编程（3）

经历了前面两个小挑战，你应该对R有点理解了。我们继续推进，今天的问题有点点复杂，复杂的不是R，而是一个数学概念：质数和质因子。任何一个合数都可以被几个质数所分解，这个性质很重要，我们将用它来解决Project Euler的第三个问题。还是和之前一样的，你需要自己在R控制台中敲打下面这些命令，根据结果自行揣摩其用处。

# 预备练习，学习for循环、建立自定义函数和其它一些函数
 
for (n in 1:10) {
    print(sqrt(n))
}
 
x <- c('hello','world','I','love','R')
for (n in x) {
    print(n)
}

笨办法学R编程（2）

本例将介绍R语言中的while循环和if条件。最终用它来解决Project Euler的第二个问题。除了练习之外你还需要了解一些斐波纳契数列的知识。废话不多说了，打开R控制台，跟着输入下面的代码，自行琢磨吧。

# 预备练习，while循环和if判断
x <- 1:10
print(x)
print(x[10])
print(x[-10])
 
i <- 1
while (i <= 10) {
  print(x[i])
  i <- i + 1
}

笨办法学R编程（1）

在倚天屠龙记中，有一人唤作火工头陀。此人练功不靠心法，只靠模仿他人招式，由外而内，自成一家。练习编程也有如此的法门，不看文字描述，只观察和模仿别人的代码。这样也可以由外而内学会编程。《笨办法学python》的作者Zed Shaw 就说过这种笨办法入门其实更简单。阳志平在他的文章《如何学习一门新的编程语言》中也讲到，初学编程要在学习区刻意的大量练习，少看理论书。

TED上一位教育家同样谈到这么一个故事，他把一个计算机扔在一个偏远的印度小村子里不去管它，在那里没有上过学的小孩就能自己学会英语和计算机的用法。实际上人脑是非常善于自我探索和学习的。因此本系列教程的特点就是只有演示代码加少量注释。通过反复模仿和练习，揣摩代码的变化和结果，你就能自行领悟其含义，并打下坚实的编程基础。

本系列每篇文章的目的都是用R语言编程来解决一个Project Euler的问题。Project Euler是一系列由易到难的计算机编程挑战，它提供了一个平台来激发我们解决问题的灵感和思路。本人写这个教程的目的有三：一是为了好玩，二是提高编程水平，三是示范说明以提供给需要的R初学者。另外从R-Blogger上了解，已经有两位高人用R在计算Project Euler，各位也可以参照他们的文章（博客1、博客2）。

Let's Go

用R来进行布丰投针实验

在3月14日也就圆周率日那一天，果壳网推出一篇文章《圆周率日特献：π究竟牛B在哪里？》。其中就提到了布丰（Buffon）用投针实验来计算π的近似值。不过这篇文章并没有详细说明如何用软件来做这个实验，也看到网上有朋友在问具体的过程。所以本文尝试用R来实现这个实验，算是狗尾续貂以作消遣。另外，Maxtrix67和魏太云都对这个问题有过精彩的论述，有兴趣的朋友可以参考一下。

简单来讲，投针实验是指假设有两根平行的线，它们之间的距离是1。随意抛掷一根长度为0.5的针，那么投针便有机会与平行线相交。如果总的投掷次数为N，发生相交的次数为X，那么可以用N/X来估计π的值。具体的故事和公式还可以参考这个文章。

R语言编程入门之七：程序查错（完）

写程序难免会出错，有时候一个微小的错误需要花很多时间来调试程序来修正它。所以掌握必要的调试方法能避免很多的无用功。

基本的除错方法是跟踪重要变量的赋值情况。在循环或条件分支代码中加入显示函数能完成这个工作。例如cat('var',var,'\n')。在确认程序运行正常后，可以将这行代码进行注释。好的编程风格也能有效的减少出错的机会。在编写代码时先写出一个功能最为简单的部分，然后在此基础上逐步添加其它复杂的功能。对输出结果进行绘图或统计汇总也能揭示一些潜在的问题。

R语言编程入门之六：循环与条件

循环
for (n in x) ｛expr}
R中最基本的是for循环，其中n为循环变量，x通常是一个序列。n在每次循环时从x中顺序取值，代入到后面的expr语句中进行运算。下面的例子即是以for循环计算30个Fibonacci数。
x <- c(1,1)
for (i in 3:30) {
    x[i] <- x[i-1]+x[i-2]
}

while (condition) {expr}
当不能确定循环次数时，我们需要用while循环语句。在condition条件为真时，执行大括号内的expr语句。下面即是以while循环来计算30个Fibonacci数。
x <- c(1,1)
i <- 3
while (i <= 30) {
    x[i] <- x[i-1]+x[i-2]
    i <- i +1
}

R语言编程入门之五：向量化运算

和matlab一样，R语言以向量为基本运算对象。也就是说，当输入的对象为向量时，对其中的每个元素分别进行处理，然后以向量的形式输出。R语言中基本上所有的数据运算均能允许向量操作。不仅如此，R还包含了许多高效的向量运算函数，这也是它不同于其它软件的一个显著特征。向量化运算的好处在于避免使用循环，使代码更为简洁、高效和易于理解。本文来对apply族函数作一个简单的归纳，以便于大家理解其中的区别所在。

所谓apply族函数包括了apply,sapply,lappy,tapply等函数，这些函数在不同的情况下能高效的完成复杂的数据处理任务，但角色定位又有所不同。

R语言编程入门之四：字符串处理

尽管R语言的主要处理对象是数字，而字符串有时候也会在数据分析中占到相当大的份量。特别是在文本数据挖掘日趋重要的背景下，在数据预处理阶段你需要熟练的操作字符串对象。当然如果你擅长其它的处理软件，比如Python，可以让它来负责前期的脏活。

获取字符串长度：nchar()能够获取字符串的长度，它也支持字符串向量操作。注意它和length()的结果是有区别的。

字符串粘合：paste()负责将若干个字符串相连结，返回成单独的字符串。其优点在于，就算有的处理对象不是字符型也能自动转为字符型。

R语言编程入门之三：输入与输出

如同ATM机一样，你首先得输入银行卡，才能输出得到钞票。数据分析也是如此，输入输出数据在分析工作中有重要的地位。下面对R语言中一些重要的输入输出函数进行小结，而其它的函数请参考官方指南。

1 读取键盘输入
如果只有很少的数据量，你可以直接用变量赋值输入数据。若要用交互方式则可以使用readline()函数输入单个数据，但要注意其默认输入格为字符型。scan()函数中如果不加参数则也可以用来手动输入数据。如果加上文件名则是从文件中读取数据。

2 读取表格文件
读取本地表格文件的主要函数是read.table()，其中的file参数设定了文件路径，注意路径中斜杠的正确用法（如"C:/data/sample.txt"），header参数设定是否带有表头。sep参数设定了列之间的间隔方式。该函数读取数据后将存为data.frame格式，而且所有的字符将被转为因子格式，如果你不想这么做需要记得将参数stringsAsFactors设为FALSE。与之类似的函数是read.csv()专门用来读取csv格式。

用R实现生命游戏(Game of Life)

生命游戏是英国数学家John Horton Conway在1970年发明的细胞自动机（cellular automaton）。它最初于1970年10月在《科学美国人》杂志中出现。生命游戏是在一个二维矩形世界中，这个世界中的每个方格居住着一个细胞。细胞的生命状态只有活着或死亡两种。细胞在下一个时刻的生死取决于相邻八个方格中活着细胞的数量。如果相邻方格活着的细胞数量过多，这个细胞会因为资源匮乏而在下一个时刻死去；相反，如果周围活细胞过少，这个细胞会因太孤单而死去。

其中一种较为丰富的设定有如下四条：

• 若相邻的存活细胞数少于2，则细胞死于寂寞
• 若相邻的存活细胞数多于3，则细胞死于拥挤
• 若相邻的存活细胞数为2或3，则细胞继续存活
• 若相邻的存活细胞数正好为3，死亡的细胞得到复活

R语言编程入门之二：对象和类

R是一种基于对象(Object)的语言，所以你在R语言中接触到的每样东西都是一个对象，一串数值向量是一个对象，一个函数是一个对象，一个图形也是一个对象。基于对象的编程(OOP)就是在定义类的基础上，创建与操作对象。

对象中包含了我们需要的数据，同时对象也具有很多属性(Attribute)。其中一种重要的属性就是它的类(Class)，R语言中最为基本的类包括了数值(numeric)、逻辑(logical)、字符(character)、列表(list)，在此基础上构成了一些复合型的类，包括矩阵(matrix)、数组(array)、因子(factor)、数据框(dataframe)。除了这些内置的类外还有很多其它的，用户还可以自定义新的类，但所有的类都是建立在这些基本的类之上的。

R语言编程入门之一：导论

简单来讲，编程是借助计算机来解决某个问题。学习编程的就是训练我们解决问题的能力。有这样一种说法：在未来，不会编程的人即是文盲。

1 为什么要学习R编程
大部分情况下解决某些问题还需要依赖一些事实或数据，结合数据分析的框架和计算工具来帮助我们决策和判断。这时候R语言编程就会派上用场。例如从大的方面来看，投资方要决定在何处建立风力发电场，就需要采集天气数据加以建模分析，评估各项目方案。从小的方面来看，个人是否应该购买某个理财产品，你需要获取过去的市场信息，模拟未来可能的变化，计算该项资产未来的期望收益和标准差。所以说学习R编程就是学习在数据环境中解决问题，从中磨练技术、锻炼智力，还能得到满足的快感。

• 学会R编程，才能理解高手的代码，并从中领会其用意并成为真正的高手。
• 学会R编程，才能深入了解函数背后的理论，从而进一步解决从未有过的新问题。

转载：来自 Google 的 R 语言编码风格指南

本文转自Xiao Nan的博客

R语言是一门主要用于统计计算和绘图的高级编程语言. 这份 R 语言编码风格指南旨在让我们的 R 代码更容易阅读、分享和检查. 以下规则系与 Google 的 R 用户群体协同设计而成.

概要: R编码风格约定

1. 文件命名: 以 .R (大写) 结尾
2. 标识符命名: variable.name, FunctionName, kConstantName
3. 单行长度: 不超过 80 个字符
4. 缩进: 两个空格, 不使用制表符
5. 空白
6. 花括号: 前括号不折行写, 后括号独占一行
7. 赋值符号: 使用 <-, 而非 =
8. 分号: 不要用
9. 总体布局和顺序
10. 注释准则: 所有注释以 # 开始, 后接一个空格; 行内注释需要在 # 前加两个空格
11. 函数的定义和调用
12. 函数文档
13. 示例函数
14. TODO 书写风格: TODO(您的用户名)