-
-
Notifications
You must be signed in to change notification settings - Fork 269
RStan Getting Started (简体中文)
RStan 是以 R 为操作界面来运行Stan的R套件。若欲获取更多有关Stan及其编程语言的信息,请移驾Stan官方网页:
以下安装指引可以适用于较早的RStan版本。在此之前,你需要 3.4.0 或更高版本的R。如有需要,您可以从这里安装更早版本的Rstan。
此外,1.2.x 或更高版本的RStudio提供了更出色的Stan支持,我们强烈推荐您使用。
假如您不确定之前是否曾经安装过Rstan,您可以首先移除旧版本:
remove.packages("rstan")
if (file.exists(".RData")) file.remove(".RData")
并且重新启动R。
在绝大多数情况下,安装Rstan仅需在R中键入
install.packages("rstan", repos = "https://cloud.r-project.org/", dependencies = TRUE)
即可。
假如您使用Linux系统,或者getOption("pkgType")
被设置为"source"
,又或者R提示您必须从source中安装最新版本的Rstan,请根据您所使用的操作系统并参考下列链接安装Rstan:
为了校验C++工具链(toolchain)是否正确安装,您可以在 RStudio或者R中, 执行
pkgbuild::has_build_tools(debug = TRUE)
如果返回 TRUE
, 则表明您的C++工具链安装正确,您可以直接跳到下一节。
否则,如果返回 FALSE
,您需要以下进一步配置:
- 若你在Windows上使用RStudio,您将会看见一个提示框提醒您安装Rtools,请点击Yes以安装。
- 若你在Windows上使用R而非RStudio,R面板中将会有类似信息提示你您安装Rtools。 请参考 如何下载安装rtools。
- 若您使用Mac, 请不要点击跳出的链接. 请点击 这里来继续安装。
- 若您使用Linux (包括Windows中的Linux子系统),请点击 这里.
如果以上步骤仍旧无法成功安装,请您移步Stan讨论区。请确保您在讨论区发帖的时候告诉我们您的操作系统,R配置环境,以及详细的输出信息。
尽管并非必须,但某些情况下优化C++配置将会提升Stan的运行速度。您仅需在初次运行时在R中复制以下命令:
dotR <- file.path(Sys.getenv("HOME"), ".R")
if (!file.exists(dotR)) dir.create(dotR)
M <- file.path(dotR, ifelse(.Platform$OS.type == "windows", "Makevars.win", "Makevars"))
if (!file.exists(M)) file.create(M)
cat("\nCXX14FLAGS=-O3 -march=native -mtune=native",
if( grepl("^darwin", R.version$os)) "CXX14FLAGS += -arch x86_64 -ftemplate-depth-256" else
if (.Platform$OS.type == "windows") "CXX11FLAGS=-O3 -march=native -mtune=native" else
"CXX14FLAGS += -fPIC",
file = M, sep = "\n", append = TRUE)
但若将optimization level设置为 O3
,您或许会在Rstan以外的R包中遭遇意外出错。在极个别的情况下,设置 -march=native -mtune=native
也将会导致Rstan出错。如果您需要改变您的 C++工具链配置,请使用
M <- file.path(Sys.getenv("HOME"), ".R", ifelse(.Platform$OS.type == "windows", "Makevars.win", "Makevars"))
file.edit(M)
如您已经依照上述安装步骤成功安装了RStan,请继续参考以下的说明内容:
RStan套件的名称为rstan (皆为小写),可以在R控制台中以library("rstan")
来进行读取
library("rstan") # 读取后可以看到相关的启动消息
您将看到启动消息弹出的提示信息。假如您所使用的电脑系统具有多核中央处理器(CPU)并且有充足的内存(RAM)来进行平行运算仿真——通常情况下没有理由不——则建议使用下列的设置步骤:
options(mc.cores = parallel::detectCores())
来并行运算运行多个马尔可夫链。
您还应键入
rstan_options(auto_write = TRUE)
俾使编译进程得自动保存。
若您使用Windows, 您或需额外配置
Sys.setenv(LOCAL_CPPFLAGS = '-march=native')
但若已上节中正确配置C++,此条并非必须。
这是在Gelman等于2003年的Bayesian Data Analysis这本书第五章第五节所提供的一个简单范例,探讨八所学校在不同的教学方式上对学术水准测验考试(SAT, Scholastic Aptitude Test)所产生的不同影响。
在此,我们将此范例简称为八校(Eight Schools)。首先编写一Stan程序的仿真模型并保存及命名为8schools.stan
在您当前的工作路径里。 若您使用RStudio 1.2.x 以上版本, 您可以点击 File -> New File -> Stan File 来新建一个stan文件。
// 存盘为 8schools.stan
data {
int<lower=0> J; // 学校数目
real y[J]; // 测试效果的预测值
real<lower=0> sigma[J]; // 测试效果的标准差
}
parameters {
real mu;
real<lower=0> tau;
real eta[J];
}
transformed parameters {
real theta[J];
for (j in 1:J)
theta[j] = mu + tau * eta[j];
}
model {
target += normal_lpdf(eta | 0, 1);
target += normal_lpdf(y | theta, sigma);
}
请确保最后一行是空白。
在此模式中,我们让theta
成为mu
及 eta
的转换参数,而不是直接的定义theta
为一模式参数。
这样的参数化过程可以更高效地运行(细节解释可参考这里)。接着我们可以输入研究数据:
schools_dat <- list(J = 8,
y = c(28, 8, -3, 7, -1, 1, 18, 12),
sigma = c(15, 10, 16, 11, 9, 11, 10, 18))
之后我们可以藉由如下R指令来得到拟合分析结果。值得注意的是,file =
必须正确指向8schools.stan
这个文件的所在位置,或者文件已经放置在目前的R工作目录中。
fit <- stan(file = '8schools.stan', data = schools_dat,
iter = 1000, chains = 4)
虽然也可使用stan
内置的model_code
功能来定义Stan模型字符串,但我们建议使用上述的方法,因为使用定义的文件在file =
这个功能可以让您输出状态描述(print statement)并反参考(dereference)其中的错误消息。
透过stan
函数所获得的输出结果fit
即为stanfit
中的S4项目 。其他方法如print
、plot
及pairs
均与拟合结果有关,故而我们可以透先前保存的fit
来检试结果。
print
提供模式参数以及对数化的后验概率密度分布(log-posterior)lp__
的估计结果总览(参考以下输出范例)。关于stanfit
的更多操作方法及细节可以进一步参考stanfit
的说明页面。
另外,我们可以使用extract
这个函数对stanfit
的项目来查看分析内容的详细结果。extract
提取了所有stanfit
项目中的抽样结果,并以条列(list)的以数组(array)方式列出参数的抽样值。
同时,S3的函数库如as.array
、as.matrix
及as.data.frame
可进一步定义stanfit
的项目为数组、矩阵及数据表(建议透过help("as.array.stanfit")
查看R的辅助文档)。
print(fit)
plot(fit)
pairs(fit, pars = c("mu", "tau", "lp__"))
la <- extract(fit, permuted = TRUE) # 回传数组列表
mu <- la$mu
### 回传三维的数组:叠代、马尔可夫链及参数
a <- extract(fit, permuted = FALSE)
### 对stanfit项目使用S3函数
a2 <- as.array(fit)
m <- as.matrix(fit)
d <- as.data.frame(fit)
鼠群也是一个著名的范例。这个例子同样也在OpenBUGS的范例中使用过,其主要是Gelfand等人于1990年所发表的研究,针对30只大鼠于五周内的成长纪录。
以下的表格表示了老鼠的代号,x
表示数据纪录时的观测天数。我们可以由链接内的数据rats.txt及Stan模式原代码rats.stan来进行范例演练。
Rat | x=8 | x=15 | x=22 | x=29 | x=36 | Rat | x=8 | x=15 | x=22 | x=29 | x=36 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | 151 | 199 | 246 | 283 | 320 | 16 | 160 | 207 | 248 | 288 | 324 | |
2 | 145 | 199 | 249 | 293 | 354 | 17 | 142 | 187 | 234 | 280 | 316 | |
3 | 147 | 214 | 263 | 312 | 328 | 18 | 156 | 203 | 243 | 283 | 317 | |
4 | 155 | 200 | 237 | 272 | 297 | 19 | 157 | 212 | 259 | 307 | 336 | |
5 | 135 | 188 | 230 | 280 | 323 | 20 | 152 | 203 | 246 | 286 | 321 | |
6 | 159 | 210 | 252 | 298 | 331 | 21 | 154 | 205 | 253 | 298 | 334 | |
7 | 141 | 189 | 231 | 275 | 305 | 22 | 139 | 190 | 225 | 267 | 302 | |
8 | 159 | 201 | 248 | 297 | 338 | 23 | 146 | 191 | 229 | 272 | 302 | |
9 | 177 | 236 | 285 | 350 | 376 | 24 | 157 | 211 | 250 | 285 | 323 | |
10 | 134 | 182 | 220 | 260 | 296 | 25 | 132 | 185 | 237 | 286 | 331 | |
11 | 160 | 208 | 261 | 313 | 352 | 26 | 160 | 207 | 257 | 303 | 345 | |
12 | 143 | 188 | 220 | 273 | 314 | 27 | 169 | 216 | 261 | 295 | 333 | |
13 | 154 | 200 | 244 | 289 | 325 | 28 | 157 | 205 | 248 | 289 | 316 | |
14 | 171 | 221 | 270 | 326 | 358 | 29 | 137 | 180 | 219 | 258 | 291 | |
15 | 163 | 216 | 242 | 281 | 312 | 30 | 153 | 200 | 244 | 286 | 324 |
y <- as.matrix(read.table('https://raw.github.com/wiki/stan-dev/rstan/rats.txt', header = TRUE))
x <- c(8, 15, 22, 29, 36)
xbar <- mean(x)
N <- nrow(y)
T <- ncol(y)
rats_fit <- stan(file = 'https://raw.githubusercontent.com/stan-dev/example-models/master/bugs_examples/vol1/rats/rats.stan')
您可以运行以下由Stan所内置的函数来演练所有的BUGS范例
model <- stan_demo()
从弹出的列表清单中选择一个范例模型来练习。在第一次使用stan_demo()
时,它会先询问您是否要下载这些范例。您需选择选第一项来将这些范例保存在rstan所安装的目录下,之后便可直接读取使用而无需重复下载。
前面所描述的model
会保存所选择的stanfit
参考范例,仿真之后就可以接着使用print
、plot
、pairs
以及extract
等函数来进行数据分析。
更多关于RStan的细节数据可以参考在rstan套件中的示范(vignette)文档。请使用help(stan)
及help("stanfit-class")
来查看stan
函数以及stanfit
这个分类中的所有介绍。
您也可参考 Stan's modeling language manual关于Stan抽样、最佳化及Stan编程语言的详细内容。
另外,Stan论坛区可以用来讨论Stan的使用、张贴范例以及关于(R)Stan的问题询问。
当您有需要时,建议能提供充足的信息如下:
- 格式正确的Stan原代码
- 数据(data)
- 必要的R原代码
- 当使用
stan
时所出现的所有错误消息,可以由verbose=TRUE
及cores=1
来显示 - 关于所使用R的相关信息,可以在R中运行
sessionInfo()
来查看
- Gelman, A., Carlin, J. B., Stern, H. S., Vehtari, A., and Rubin, D. B. (2013). Bayesian Data Analysis, CRC, 3nd Edition.
- Stan Development Team. Stan Modeling Language User's Guide and Reference Manual.
- Gelfand, A. E., Hills S. E., Racine-Poon, A., and Smith A. F. M. (1990). "Illustration of Bayesian Inference in Normal Data Models Using Gibbs Sampling", Journal of the American Statistical Association, 85, 972-985.
- Stan
- R
- BUGS
- OpenBUGS
- JAGS
- Rcpp