This commit is contained in:
LSGOMYP
2021-10-13 21:22:16 +08:00
22 changed files with 494 additions and 1 deletions

View File

@@ -0,0 +1,132 @@
# ARMA 时间序列模型与预测
### 白噪声
可以通过 Box-Ljung 检验来检验序列是否为白噪声:
```R
> set.seed(100)
> data = rnorm(100)
> Box.test(data, type='Ljung', lag = log(length(data)))
Box-Ljung test
data: data
X-squared = 8.896, df = 4.6052, p-value = 0.09169
```
从结果中可以看见 $p = 0.09169 > 0.05$,因此无法拒绝序列为白噪声的假设。下面绘制一下该序列的图像以及 ACF 图像:
```R
> op <- par(mfrow=c(2, 1), mar=c(5, 4, 2, 2) + .1)
> plot(ts(data))
> acf(z, main = "")
> par(op)
```
![image-20210924205843412](.\typora-user-images\image-20210924205843412.png)
### AR(p) 序列
先给出一个 AR(1) 的手工计算例子,在这个例子中可以通过 $y_{n-1}$ 预测 $y_n$
```R
> n <- 200
> x <- rnorm(n)
> f = c(1,2)
> y1 = x[1]; y1
[1] 1
> y2 = x[2] + f[1] * y1; y2
[1] 3
> y3 = x[3] + f[1] * y2 + f[2] * y1; y3
[1] 8
> y4 = x[4] + f[1] * y3 + f[2] * y2; y4
[1] 18
```
实际上也可以通过 `filter()` 函数来计算:
```R
> x = 1:10; x
[1] 1 2 3 4 5 6 7 8 9 10
> y = filter(x, c(1, 2), method='r'); y
Time Series:
Start = 1
End = 10
Frequency = 1
[1] 1 3 8 18 39 81 166 336 677 1359
```
下面计算一个 $p=3$ 的例子:
```R
> n <- 200
> x <- rnorm(n)
> f = c(.3, -.7, .5)
> y <- rep(0, n)
> y[1:3] = x[1:3]
> for (i in 4:n) {
+ y[i] <- f[1]*y[i-1] +f[2]*y[i-2] + f[3]*y[i-3] + x[i]
+ }
> op <- par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)
> plot(ts(y), xlab="", ylab="AR(3)")
> acf(y, main="", xlab="")
> pacf(y, main="", xlab="")
> par(op)
```
![image-20210923195745948](.\typora-user-images\image-20210923195745948.png)
同样地,也可以通过 `filter()` 函数来计算:
```R
> y = filter(x, c(.3, -.7, .5), method='r'); y
> op <- par(mfrow=c(3,1), mar=c(2,4,2,2)+.1)
> plot(ts(y), xlab="", ylab="AR(3)")
> acf(y, main="", xlab="")
> pacf(y, main="", xlab="")
```
![image-20210923200033331](.\typora-user-images\image-20210923200033331.png)
结果是一样的。
### MA 模型
直接使用 `filter()` 函数计算 MA 模型:
```R
> x = 1:10
> y = filter(x, filter = c(.5, .3)); y # filter 函数未设置 method默认为'c',即使用滑动平均
Time Series:
Start = 1
End = 10
Frequency = 1
[1] 1.3 2.1 2.9 3.7 4.5 5.3 6.1 6.9 7.7 NA
```
作业:编写手工计算 MA 模型的代码
### ARIMA 模型相关
R 语言中自带的 `arima.sim()` 函数可以模拟生成 AR、MA、ARMA 或 ARIMA 模型的数据。其原型为:
```R
arima.sim(model, n, rand.gen = rnorm, innov = rand.gen(n, ...),
n.start = NA, start.innov = rand.gen(n.start, ...),
...)
```
其中,`model` 是一个列表,用于指定各模型的系数;`order` 是 ARIMA(p, d, q) 中 $(p, d, q)$ 三个元素的向量,$p$ 为 AR 阶数, $q$ 是 MA 的阶数,$d$ 是差分阶数。例如,模拟如下的 ARIMA(1, 1, 1) 模型并产生长度为300的样本
$Y_t = X_t - X_{t-1}, X_t = -0.9X_{t-1} + \varepsilon_t + 0.5\varepsilon_{t-1}, \varepsilon_t \sim WN(0, 2^2)$
R 语言代码为:
```R
> x <- 2.0 * arima.sim(model = list(
+ ar = c(-0.9), ma = c(0.5), order = c(1, 1, 1)), n=300)
```

View File

@@ -0,0 +1,22 @@
<EFBFBD><EFBFBD><EFBFBD><EFBFBD>,ָ<EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD>Ȼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>,ɽ<EFBFBD><EFBFBD><EFBFBD>߿<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>,ɽ<EFBFBD><EFBFBD>ʡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>,<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>˿<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
1,1999<EFBFBD><EFBFBD>,3203.6,51.03,15.93<EFBFBD><EFBFBD>,19.45,6.07<EFBFBD><EFBFBD>,9.86<EFBFBD><EFBFBD>,,,
2,2000<EFBFBD><EFBFBD>,3247.8,47.72,13.25<EFBFBD><EFBFBD>,18.59,5.77<EFBFBD><EFBFBD>,7.48<EFBFBD><EFBFBD>,,,
3,2001<EFBFBD><EFBFBD>,3271.6,42.57,13.06<EFBFBD><EFBFBD>,19.24,5.90<EFBFBD><EFBFBD>,7.16<EFBFBD><EFBFBD>,,,
4,2002<EFBFBD><EFBFBD>,3293.7,42.22,12.86<EFBFBD><EFBFBD>,20.16,6.14<EFBFBD><EFBFBD>,6.82<EFBFBD><EFBFBD>,,,
5,2003<EFBFBD><EFBFBD>,3314.3,40.52,12.26<EFBFBD><EFBFBD>,19.94,6.04<EFBFBD><EFBFBD>,6.22<EFBFBD><EFBFBD>,,,
6,2004<EFBFBD><EFBFBD>,3335.1,41.1,12.36<EFBFBD><EFBFBD>,20.32,6.11<EFBFBD><EFBFBD>,6.25<EFBFBD><EFBFBD>,,,
7,2005<EFBFBD><EFBFBD>,3355.2,40.21,12.02<EFBFBD><EFBFBD>,20.07,6.00<EFBFBD><EFBFBD>,6.02<EFBFBD><EFBFBD>,,,
8,2006<EFBFBD><EFBFBD>,3374.6,38.63,11.48<EFBFBD><EFBFBD>,19.28,5.73<EFBFBD><EFBFBD>,5.75<EFBFBD><EFBFBD>,,,
9,2007<EFBFBD><EFBFBD>,3392.6,38.26,11.31<EFBFBD><EFBFBD>,20.23,5.98<EFBFBD><EFBFBD>,5.33<EFBFBD><EFBFBD>,,,
10,2008<EFBFBD><EFBFBD>,3410.6,38.51,11.32<EFBFBD><EFBFBD>,20.44,6.01<EFBFBD><EFBFBD>,5.31<EFBFBD><EFBFBD>,37,18.43,18.57
11,2009<EFBFBD><EFBFBD>,3427.4,37.16,10.87<EFBFBD><EFBFBD>,20.45,5.98<EFBFBD><EFBFBD>,4.89<EFBFBD><EFBFBD>,37.2,17.63,19.57
12,2010<EFBFBD><EFBFBD>,3574.1,38.06,10.68<EFBFBD><EFBFBD>,19.18,5.38<EFBFBD><EFBFBD>,5.30<EFBFBD><EFBFBD>,36.4,18.44,17.96
13,2011<EFBFBD><EFBFBD>,3593.3,37.5,10.47<EFBFBD><EFBFBD>,20.1,5.61<EFBFBD><EFBFBD>,4.86<EFBFBD><EFBFBD>,34.2,18.46,15.74
14,2012<EFBFBD><EFBFBD>,3610.8,38.53,10.70<EFBFBD><EFBFBD>,20.99,5.83<EFBFBD><EFBFBD>,4.87<EFBFBD><EFBFBD>,33.9,20.81,13.09
15,2013<EFBFBD><EFBFBD>,3629.8,39.15,10.81<EFBFBD><EFBFBD>,20.17,5.57<EFBFBD><EFBFBD>,5.24<EFBFBD><EFBFBD>,36.1,21.53,14.57
16,2014<EFBFBD><EFBFBD>,3648,39.74,10.92<EFBFBD><EFBFBD>,21.58,5.93<EFBFBD><EFBFBD>,4.99<EFBFBD><EFBFBD>,34.2,21.44,12.76
17,2015<EFBFBD><EFBFBD>,3664.1,36.49,9.98<EFBFBD><EFBFBD>,20.33,5.56<EFBFBD><EFBFBD>,4.42<EFBFBD><EFBFBD>,34.2,22.27,11.93
18,2016<EFBFBD><EFBFBD>,3681.6,37.79,10.29<EFBFBD><EFBFBD>,20.27,5.52<EFBFBD><EFBFBD>,4.77<EFBFBD><EFBFBD>,33.9,22,11.9
19,2017<EFBFBD><EFBFBD>,3702.4,40.83,11.06<EFBFBD><EFBFBD>,20.12,5.45<EFBFBD><EFBFBD>,5.61<EFBFBD><EFBFBD>,31.7,22.14,9.56
20,2018<EFBFBD><EFBFBD>,3718.3,35.73,9.63<EFBFBD><EFBFBD>,19.74,5.32<EFBFBD><EFBFBD>,4.31<EFBFBD><EFBFBD>,30.5,22.39,8.11
21,2019<EFBFBD><EFBFBD>,3729.2,33.97,9.12<EFBFBD><EFBFBD>,21.79,5.85<EFBFBD><EFBFBD>,3.27<EFBFBD><EFBFBD>,31.4,25.34,6.06
1 序号 指标 总人数 出生人数 出生率 死亡人数 死亡率 自然增长率 山西高考人数 山西省招生人数 教育人口流出
2 1 1999年 3203.6 51.03 15.93‰ 19.45 6.07‰ 9.86‰
3 2 2000年 3247.8 47.72 13.25‰ 18.59 5.77‰ 7.48‰
4 3 2001年 3271.6 42.57 13.06‰ 19.24 5.90‰ 7.16‰
5 4 2002年 3293.7 42.22 12.86‰ 20.16 6.14‰ 6.82‰
6 5 2003年 3314.3 40.52 12.26‰ 19.94 6.04‰ 6.22‰
7 6 2004年 3335.1 41.1 12.36‰ 20.32 6.11‰ 6.25‰
8 7 2005年 3355.2 40.21 12.02‰ 20.07 6.00‰ 6.02‰
9 8 2006年 3374.6 38.63 11.48‰ 19.28 5.73‰ 5.75‰
10 9 2007年 3392.6 38.26 11.31‰ 20.23 5.98‰ 5.33‰
11 10 2008年 3410.6 38.51 11.32‰ 20.44 6.01‰ 5.31‰ 37 18.43 18.57
12 11 2009年 3427.4 37.16 10.87‰ 20.45 5.98‰ 4.89‰ 37.2 17.63 19.57
13 12 2010年 3574.1 38.06 10.68‰ 19.18 5.38‰ 5.30‰ 36.4 18.44 17.96
14 13 2011年 3593.3 37.5 10.47‰ 20.1 5.61‰ 4.86‰ 34.2 18.46 15.74
15 14 2012年 3610.8 38.53 10.70‰ 20.99 5.83‰ 4.87‰ 33.9 20.81 13.09
16 15 2013年 3629.8 39.15 10.81‰ 20.17 5.57‰ 5.24‰ 36.1 21.53 14.57
17 16 2014年 3648 39.74 10.92‰ 21.58 5.93‰ 4.99‰ 34.2 21.44 12.76
18 17 2015年 3664.1 36.49 9.98‰ 20.33 5.56‰ 4.42‰ 34.2 22.27 11.93
19 18 2016年 3681.6 37.79 10.29‰ 20.27 5.52‰ 4.77‰ 33.9 22 11.9
20 19 2017年 3702.4 40.83 11.06‰ 20.12 5.45‰ 5.61‰ 31.7 22.14 9.56
21 20 2018年 3718.3 35.73 9.63‰ 19.74 5.32‰ 4.31‰ 30.5 22.39 8.11
22 21 2019年 3729.2 33.97 9.12‰ 21.79 5.85‰ 3.27‰ 31.4 25.34 6.06

View File

@@ -0,0 +1,3 @@
#### 本部分存储了ARMA模型案例山西省教育人口逆差数据与其他相关数据可以用stata或R语言进行分析。
#### 数据来源:人口普查数据源于国家统计局;山西省人口数据来源山西省统计局《统计年鉴》。

View File

@@ -0,0 +1,120 @@
# R 语言基础
R语言是一门常用于数据分析、统计建模的计算机语言它与主流的C/C++、Java、Python等语言相比支持更多的数据类型例如向量、矩阵同时提供了多种统计和数学计算方法。
可以前往 https://www.r-project.org/ 下载R语言解释器并且推荐使用 RStudio 这个R语言的集成开发环境。RStudio 可以在 https://www.rstudio.com/ 下载。
## 四则运算
R语言的四则运算与大部分语言相同使用 `+`, `-`, `*`, `/`, `^` 来表示加、减、乘、除和乘方。数值可以写成 `123`, `-123`, `123.45`, `1.23E-5`这样的形式。其中`1.23E-5`表示$1.23 \times 10^{-5}$ 。
## 字符串
用单引号 `'` 或双引号`"` 包裹起来的文字内容为字符串,如 `'hello, world'``'123456'`
## 向量
R 语言中可以通过 `c(...)` 来声明一个向量,例如 $\boldsymbol{x} = (1, 2, 3, 4, 5)$ 可以通过 `x <- c(1, 2, 3, 4, 5)` 来声明。
对两个长度相等的向量进行四则运算的效果是向量中的每一个元素都与另一个向量中的每一个元素进行四则运算。而一个向量与一个数进行四则运算的效果是该向量中的所有元素都与这个数进行四则运算。如:
```R
> x <- c(1, 2, 3, 4, 5)
> y <- c(1, 2, 3, 4, 5)
> x + y # 输出: [1] 2 4 6 8 10
> x + 1 # 输出: [1] 2 3 4 5 6
```
## 矩阵
R 语言中可以通过 `matrix()` 函数来创建矩阵。`matrix()` 的原型为 `matrix(data=NA, nrow=1, ncol=1, byrow=FALSE, dinames=NULL)` 其中:
- `data` : 矩阵的元素,通常为向量。
- `nrow``ncol` : 设定矩阵的行数和列数,一般只需设定其一,另一个会根据数据长度算出。
- `byrow` : 设定矩阵的填充方式,值为 `TRUE` 时按行填充。默认为 `FALSE` ,即按列填充。
通过 `matrix()` 创建矩阵的例子如下:
```R
> mat1 <- matrix(1:6, nrow=2) # 元素为1到6 两行,按列填充
> mat1
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> mat2 <- matrix(1:6, nrow=2, byrow=TRUE)
> mat2
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
```
此外还可以使用 `dim()` 通过向量来创建矩阵,例如:
```R
> mat3 <- 1:6
> dim(mat3) <- c(3, 2) # 将元素为1到6变为3行2列的矩阵
> mat3
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
```
## 生成时间序列
通过 `ts()` 函数可以将一个向量或矩阵转成一个一元或多元的时间序列ts对象`ts()` 函数的原型为:`ts(data = NA, start = 1, end = numeric(0), frequency = 1, deltat = 1, ts.eps = getOption("ts.eps"), class, names)`。其中:
- `data` 要生成时间序列的向量或矩阵。
- `start` 第一个观测值的时间。
- `end` 最后一个观测值的时间。
- `frequency` 单位时间内观测值的频数。
- `deltat` 两个观测值之间的时间间隔。
- `ts.eps` 序列之间的误差限。若序列之间的频率差异小于 `ts.eps` 则认为这些序列的频率相等。
- `class` 对象的类型。一元序列默认为 `ts`,多元序列默认为 `c("mts", "ts")`
- `names` 给出多元序列中每个一元序列的名称,默认为 `Series 1, Series 2, ...`
下面是一个例子:
```R
> ts(1:26, start = 1986)
Time Series:
Start = 1986
End = 2011
Frequency = 1
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
[18] 18 19 20 21 22 23 24 25 26
```
## 画图
R 语言中可以直接使用 `plot()` 函数来绘制图像,例如:
```R
> x <- rnorm(10)
> plot(x)
> x
[1] -0.3329234 1.3631137 -0.4691473 0.8428756
[5] -1.4579937 -0.4003059 -0.7764173 -0.3692965
[9] 1.2401015 -0.1074338
```
![image-20210924211135207](.\typora-user-images\image-20210924211135207.png)
还可以直接绘制时间序列的图像:
```R
> plot(ts(x))
```
![image-20210924211220048](.\typora-user-images\image-20210924211220048.png)
为了将多幅图画在一起,可以使用 `par()` 函数:
```R
> op <- par(mfrow=c(2, 1), mar=c(5, 4, 2, 2) + .1) # mfrow 指定了图像矩阵为两行一列即画两幅图每行一幅mar 指定了图像的边界,分别是下、左、上、右,可以根据喜好指定
> plot(ts(x))
> acf(z, main = "") # 计算 acf 函数
```
![image-20210924211736213](.\typora-user-images\image-20210924211736213.png)

Binary file not shown.

View File

@@ -1,4 +1,34 @@

## 简介
岳昆,你们自己来填写哈
时间序列也称动态序列,是指将某种现象的指标数值按照时间顺序排列而成的数值序列。时间序列数据本质上反映的是某个或某些随机变量时间不断变化的趋势,而时间序列方法的核心就是从数据中挖掘出这种规律,并利用其对将来的数据做出估计。随着社会的进步和计算机技术的发展,时间序列分析的应用越来越广泛, 已在经济、气象、地质、水文、军事等领域产生了显著的经济效益和社会效益
时间序列应用方向可大致分成三个部分,分别是描述过去、分析规律和预测未来。
### Task00前置内容2天
- 熟悉组队学习规则
- 安装 R 解释器和 RStudio
- 学习 R 语言的基础语法和常用函数
- 数理统计基础(教程第一章)
### Task01手算时间序列2天
- 学习几种移动平均和指数平滑的时间序列处理方法(教程第二章)
- 根据所学方法的原理编写相应的 R 语言代码
### Task02常用的时间序列模型3天
- 学习常用的时间序列模型(教程第三章)
### Task03ARMA 模型与预测4天
- 了解时间序列分析的流程
- 了解几种常见的时间序列及其定阶
参考资料:
- 《数学建模算法与程序》司守奎-海军航空工程学院
- 数学模型 .高等教育出版社[引用日期2019-04-14]

Binary file not shown.

After

Width:  |  Height:  |  Size: 7.7 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 7.9 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 6.7 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 7.2 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 8.8 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 8.7 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 8.2 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 8.5 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 8.5 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 5.7 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 6.7 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 8.2 KiB

View File

@@ -0,0 +1,66 @@
# 手算时间序列
## 简单移动平均
我们先来自己实现一个计算简单移动平均的函数:
```R
mySMA <- function (x, n) {
sma <- c()
sma[1:(n-1)] <- NA
for (i in n:length(x)) {
sma[i] <- mean(x[(i-n+1):i])
}
return(sma)
}
> x = c(2, 3, 3, 4, 2, 3, 3, 5, 2, 3)
> mySMA(x, n = 4)
[1] NA NA NA 3.00 3.00 3.00 3.00 3.25 3.25 3.25
```
R 语言中的 TTR 包提供的 SMA 函数也可以实现简单移动平均的计算。首先确保 R 语言解释器的版本在 4.2.0以上(在 RStudio 中键入 `R.version` 中可查看 R 解释器的版本)。进入 https://cran.r-project.org/web/packages/TTR/index.html ,下载系统对应的 xts, zoo, curl 和 TTR 包所需要的包,并在 RStudio -> Tools -> Installs Packages 中选择下载的包进行安装。随后通过 `library(xts)` , `library(zoo)`, `library(curl)``library(TTR)` 导入相应的包。
例如对于数据 $2, 3, 3, 4, 2, 3, 3, 5, 2, 3$ 。 假设当前观测值只受到过去4期数值影响可以通过如下的方法计算
```R
> x = c(2, 3, 3, 4, 2, 3, 3, 5, 2, 3)
> SMA(x, n = 4)
[1] NA NA NA 3.00 3.00 3.00 3.00 3.25 3.25 3.25
```
## 指数移动平均
老样子,我们先来试着自己实现一下指数移动平均:
```R
myEMA <- function (price,n){
ema <- c()
ema[1:(n-1)] <- NA
ema[n]<- mean(price[1:n])
beta <- 2/(n+1)
for (i in (n+1):length(price)){
ema[i]<-beta * price[i] +
(1-beta) * ema[i-1]
}
return(ema)
}
```
```R
> x = c(2, 3, 3, 4, 2, 3, 3, 5, 2, 3)
> myEMA(x, n = 4)
[1] NA NA NA 3.000000 2.600000
[6] 2.760000 2.856000 3.713600 3.028160 3.016896
```
再来试试 TTR 包提供的指数移动平均:
```R
> EMA(x, n=4)
[1] NA NA NA 3.000000 2.600000
[6] 2.760000 2.856000 3.713600 3.028160 3.016896
```

View File

@@ -0,0 +1,9 @@
# 时间序列参考网课与教材
本部分推荐的网课与教材均为辅助学习使用,不做强制打卡要求。大家可以根据实际情况来自行选择。
- 《时间序列分析及应用R语言》系列课程视频
https://www.bilibili.com/video/BV1Hz4y167fD?from=search&seid=8685454508264223591&spm_id_from=333.337.0.0
- 时间序列分析的基本思路与步骤
https://www.bilibili.com/video/BV1454y1v79B?from=search&seid=8685454508264223591&spm_id_from=333.337.0.0
- 高级计量经济学及Stata应用第二版陈 强 编著-高等教育出版社

View File

@@ -0,0 +1,111 @@
# 统计与时间序列分析基础
## 平均值
```R
> y = array(1:20, dim=c(4, 5))
> y
[,1] [,2] [,3] [,4] [,5]
[1,] 1 5 9 13 17
[2,] 2 6 10 14 18
[3,] 3 7 11 15 19
[4,] 4 8 12 16 20
> mean(y) # 对矩阵中的所有元素求平均值
[1] 10.5
> colMeans(y) # 计算每行的均值
[1] 2.5 6.5 10.5 14.5 18.5
> rowMeans(y) # 列均值
[1] 9 10 11 12
```
## 方差和标准差
```R
> x = exp(seq(-1, 3, by=0.1)) # 从-1到3步长为0.1
> x
[1] 0.3678794 0.4065697 0.4493290 0.4965853
[5] 0.5488116 0.6065307 0.6703200 0.7408182
[9] 0.8187308 0.9048374 1.0000000 1.1051709
[13] 1.2214028 1.3498588 1.4918247 1.6487213
[17] 1.8221188 2.0137527 2.2255409 2.4596031
[21] 2.7182818 3.0041660 3.3201169 3.6692967
[25] 4.0552000 4.4816891 4.9530324 5.4739474
[29] 6.0496475 6.6858944 7.3890561 8.1661699
[33] 9.0250135 9.9741825 11.0231764 12.1824940
[37] 13.4637380 14.8797317 16.4446468 18.1741454
> var(x) # 计算方差
[1] 29.35325
> sd(x) # 计算标准差
[1] 5.417864
```
## 偏度和峰度
可以通过自己编写函数来计算偏度和峰度:
```R
> skewness <- function(x) { sum(((x-mean(x)) ^ 3)) / length(x) } # 计算偏度
> skewnexx(x)
[1] 197.8397
> kurtosis <- function(x) {
+ a = mean(x)
+ n = length(x)
+ m4 = sum((x-a)^4)/n
+ m2 = sum((x-a)^2)/n
+ kurt = m4/m2^2-3
+ kurt
+ } # 计算峰度
> kurtosis(x)
[1] 0.6260693
```
也可以通过 moments 包中的 skewness 和 kurtosis 函数来计算偏度和峰度。
首先通过 `> install.packages(moments)` 来安装 moments 包,然后通过 `> library(moments) ` 来使用这个包。之后就可以直接通过 `skewness()``kurtosis()` 来计算偏度和峰度。
## 重要的概率分布
### 正态分布
通过 `rnorm()` 函数可以产生服从正态分布的随机数,`rnorm()` 函数的原型为 `rnorm(n, mean=0, sd=1)` ,其中 `n` 为生成数据个数,`mean` 为生成数据的均值,`sd` 为生成数据的标准差。例如:
```R
> norm = rnorm(1000, 0, 1)
> mean(norm)
[1] 0.001165962
> sd(norm)
[1] 1.000664
> hist(norm, prob=TRUE, 30) # 绘制频率分布直方图
> curve(dnorm, add=TRUE) # 为直方图添加正态分布曲线
```
![image-20210922175847365](.\typora-user-images\image-20210922175847365.png)
### $\chi^2$ 分布
可以通过 `rchisq(n, df)` 来产生 n 个 服从自由度为 df 的 $\chi^2$ 分布的随机数
```R
> rchisq(n=1000, df=20)
> hist(chi, prob=TRUE, 30)
```
![image-20210922180552465](.\typora-user-images\image-20210922180552465.png)
### $t$ 分布
```R
> t = rt(n=1000, df=5)
> hist(t, prob=TRUE, 30)
```
![image-20210922185944045](.\typora-user-images\image-20210922185944045.png)
### $F$ 分布
```R
> f = rf(n=1000, df1=10, df2=50)
> hist(f, prob=TRUE, 30)
```
![image-20210922191248078](.\typora-user-images\image-20210922191248078.png)