程序代写代做代考 Excel clock html graph JDBC C 有哪些常见的统计软件

有哪些常见的统计软件
它们的特点是什么

有哪些统计软件
出处: https://en.wikipedia.org/wiki /Comparison_of_statistical_pa ckages

非统计软件:Excel+VBA
• Excel不需要多介绍了吧?
• 绝大多数同学毕业后最常用的软件
• 市场占有率最大的桌面级数据制表/分 析软件
• 绝大多数统计软件都提供了Excel数据 接口
• 复旦师生可以免费下载使用
• Why VBA?
• 数据操作比较复杂 • 重复性的数据操作 • 数据的量变大
• 不是专门的统计分析软件
• 专业的事情留给专业的软件去做

统计软件:SAS
• http://www.sas.com
• 企业级产品，体积庞大，价格贵，
成熟方法几乎无所不包
• 语言比较陈旧，少量菜单操作，
但功能强，计算效率高
• 最近这些年针对新方法新应用提
供了很多更新的产品
• 业界标准，少数几家FDA认可的分 析软件
• 要充分利用其功能，一定要学会 宏语言。

统计软件:SPSS
• https://www.ibm.com/analyti cs/spss-statistics-software
• 功能全，界面友好，菜单操作 • 商业软件，复旦师生可以免费
下载使用
• Windows/Mac/Linux各个操作系 统平台都有

统计软件:Stata
• 功能全，特别是计量经济相关的模型
• 界面友好，可菜单操作，可编程计算
• 计量经济学界使用较多，但是业界使用 很少(因为两者关心的问题很多不一样)
• 商业软件，贵(最低配置学生版198美 元)
• 用户支持不友好(教材和参考书太少)
• 各个操作系统平台都有

统计软件:MATLAB
• www.mathworks.com
• 命令行操作为主
• 数值计算功能全而且强大
• 灵活，跟最新的方法跟得紧
• 各个操作系统平台都有
• 复旦师生可以免费下载使用

统计软件:Python
• 粘合语言，主要靠各种库/包实现 相关功能
• 计算机界流行
• 现在处于2.x到3.x的转移过程中，
两者之间有部分不兼容情况 • 使用人员多，用户支持好

统计软件:R
• 针对数据分析而编写的软件，其 前身S语言，曾经获得ACM大奖
• 免费软件
• 各个操作系统版本都有
• 微软公司跟进支持，有改进版本 Microsoft R Open (MRO)
• 统计学界标配
• 建议搭配Rstudio编程环境一起使 用。

统计软件:RStudio
• R的IDE(Integrated Development Environment)

统计软件在做什么

使用统计软件分析数据的流程
• 输入
• 手工输入
• 从文件输入(文本文件、二进制文件) • 从数据库输入
• 数据整理(数据清理、转换、等等)
• 分析
• 使用软件内置的分析工具直接分析
• 使用软件自带的编程语言编写程序来分析 • 与其他软件或语言结合共同分析
• 输出
• 图形或数值、文字
• 屏幕或文件
• 指定格式的数据集(可供其他软件使用)

我该如何选择统计软件
如何才能学好统计软件
该学到哪个程度

如何选择统计软件
• 通用语言与专用语言
• 不存在所谓“最好的”，只有最“合适”的
• 尺有所短寸有所长
• 如何选择?
• 行业内(课题组内)主流是什么软件?
• 使用软件的目的是什么?(应用?研发?)
• 菜单操作 v.s. 命令行操作?
• 批处理 v.s. 交互操作?
• 为什么选用R? • 版权
• 就业市场的需求
• 学生的背景知识
• 要求掌握的内容及程度

如何学好统计软件
• 自己多动手编程序
• 多读别人的程序，“拿来主义” • 学习命名风格
• 学习代码结构风格
• 学习代码背后的思维
• 学了立即就用，多用，熟能生巧
• 各种论坛请教高手
• 统计之都(COS)论坛 • CSDN网站
• StackOverFlow网站
• 注意:某些部分内容会比较琐碎，要注意及时总结(如命令的选项)
• 进阶学习内容 • 数据结构
• 算法

我该学到哪个程度
• Q1:你要用到哪个程度?
• Q2:跟你的专业有关系，跟你从事的工作有关，也跟你要跟的导师有关系
硕博(统计，管
科，信管)
硕博(非统计管 科信管)
金融硕
MPAcc
国际商务硕士
数据输入输出
手工输入
√
√
√
√
√
文件输入输出
√
√
√
√
√
数据库输入输出
√
√
数据整理
√
√
√
√
√
数据分析
内置分析工具
√
√
√
√
√
自带的编程语言
√
√
√
与其他软件或语
言结合共同分析
√
输出
图形数值文字
√
√
√
√
√
屏幕，文件
√
√
√
√
√
指定格式数据集
√
√
√
Excel
Excel
其他要学的内容
数据结构与算法

几种错误的认识和倾向
• 错误1:要学就学最牛的那个
• 解释1:不存在所谓“最好的”，只有最“合适”的
• 错误2:学了某个软件之后就可以一劳永逸
• 解释2:软件永远在变，学习永远在路上，但基本是相通的，可以
举一反三
• 错误3:习惯了某个软件就不愿意改
• 解释3:
• 人都有惰性，对某一个工具掌握了之后就会有依赖性，凡事总想拿现成的 工具去解决，而且还会找各种理由为自己辩护，说这个软件怎么怎么好， 别的软件怎么怎么差。真实的理由往往就是一个字，懒!
• 跳出舒适区，逼自己一把，否则就被别人干掉

使用统计软件分析数据的流程

使用统计软件分析数据的流程
• 实际工作中分析数据的顺序 • 输入
• 数据整理(数据清理、转换、等等) • 分析
• 输出

使用统计软件分析数据的流程
• 我们看一个小例子
• 我们要计算体重的均值，标准差，体重与年龄的相关系数，并 作出体重与年龄的散点图

使用统计软件分析数据的流程
• 实际工作中分析数据的顺序 • 输入
• 数据整理(数据清理、转换、等等) • 分析
• 输出
age = c(1, 3, 5, 2, 11, 9, 3, 9, 12, 3)
weight = c(4.4, 5.3, 7.2, 5.2, 8.5, 7.3, 6.0, 10.4, 10.2, 6.1)
mean(weight)
sd(weight)
cor(age, weight)
plot(age, weight)
q()

使用统计软件分析数据的流程
• 实际工作中分析数据的顺序 • 输入
• 数据整理(数据清理、转换、等等) • 分析
• 输出
• 讲解的顺序会有些不一样 • 输入
• 输出
• 数据操作、清理(数据清理、转换、等等) • 分析

认识数据

什么是数据?
• 信息的载体均可视为数据 • 数字
• 文本 • 图像 • 音频 • 视频 • ……
• 本课程主要关心的是数字类型的数据
• 这是最常见的数据，也是被研究得最为透彻的
• 对于其他类型的分析一般也是先被转化为数字类型的数据，再被研究和分析的。
• 数据一般长什么样?

数据的组织与表示
表头
身份证号
姓名
性别
成绩
绩点
分数
用时
NA
葛老
男
A
4
98
2
120230xxxx
闻染
女
A
4
100
2
310109xxxx
李泌
男
A
4
100
2
310108xxxx
张小敬
男
A
4
99
2
100005xxxx
龙波
男
B
3
70
2
374983xxxx
崔器
男
A
4
98
2
546284xxxx
杨玉环
女
A+
4.3
110
1.2
201251xxxx
徐宾
男
–
–
–
–
201214xxxx
王韫秀
女
B-
2.7
60
4
654521xxxx
姚汝能
男
A+
4.3
110
1
观 测 ( 个 体 )
观 测 值
变量

数据的类型与信息:定性变量(Categorial)
• 例子:性别、颜色、身份证号前6位
• 内涵信息:
• 每一个值表示一种属性
• 不同的值表示不同的属性
• 不同的值仅仅表示属性的不同，没有其他诸如先后、高低等含义
• 计算:
• 仅可以做是否相同，是否不同的逻辑运算。
• 注:
• 计算机中，常将各个值进行编码存储与计算 • 又称为名义变量(Norminal data)

数据的类型与信息:定序变量(Ordinal)
• 例子:考试成绩等级，军衔，名次
• 内涵信息:
• 每一个值表示一种属性
• 不同的值表示不同的属性
• 不同的值之间有先后、高低等含义
• 计算:
• 可以做是否相同，是否不同的逻辑运算。 • 也可以比较大小(先后)
• 注:
• 计算机中，常将各个值进行编码存储与计算

数据的类型与信息:定距变量(Interval)
• 例子:时间
• 内涵信息:
• 每一个值表示一种属性
• 不同的值表示不同的属性
• 不同的值之间有先后、高低等含义 • 不同的值之间可以计算间隔的大小
• 计算:
• 可以做是否相同，是否不同的逻辑运算。 • 也可以比较大小(先后)
• 可以做加减法

数据的类型与信息:定比变量(Proportional)
• 例子:距离，体重，时间间隔(考试用时)
• 内涵信息:
• 每一个值表示一种属性
• 不同的值表示不同的属性
• 不同的值之间有先后、高低等含义 • 不同的值之间可以计算间隔的大小 • 有绝对0点
• 计算:
• 可以做是否相同，是否不同的逻辑运算。 • 也可以比较大小(先后)
• 可以做四则运算
• 某些情况下可以做初等函数的运算

数据的类型与信息:小结
• 名义:男女、颜色
• 次序:名次、等级(军衔) • 间隔:时间
• 比率:距离、体重
=,≠
>, < 四则运算 其他运算(如指数、开方等 名义 Y N N N 次序 Y Y N N 间隔 Y Y +、- N(除非可以确定零点) 比率 Y Y Y Y ) 数据的类型(其他标准) • 按时间 • 横截面数据 • 时间序列数据 • 面板数据 • 按来源 • 现有数据(历史数据，如统计年鉴、气象记录) • 调查数据(新iphone的市场前景) • 按总体与部分 • 总体数据(描述统计) • 样本数据(描述统计 + 推断统计) • 按变化方式 • 常量:光速、e、𝜋 • 变量:日期、年龄 • 随机变量: • 复旦大学管理学院2018级本科生内部能成几对? • 你的下一个孩子是男孩还是女孩? • 今天会接到几个电话? • 这个月的网络流量是多少? 对象 什么是对象(object) • 当数据输入之后，它不仅仅只是一个数据，同时要考虑对于这 个数据集我们可以进行哪些分析和操作。 • 在R中，所有的数据与可以施加在此类数据上的属性、函数、 方法等，是一个有机的整体，我们统称为“对象”。 • 在R中，一切皆对象(object)。 对象的命名规则 • 命名规则: • 以字符开头 • 可以由字母、下划线“_”、点“.”、数字等组成 • 大小写敏感(例:A与a不能通用) • 常量名pi，特殊字符NA，NaN 等已被系统保留 • 为避免混淆，尽量避免使用系统已经保留的变量名，函数名，等等。 对象的操作命令 x=1 e1 = rnorm(100) e2 = matrix(1:12, 3, 4) ae = x + 3 m1 = 1:5 xm = "xm" objects() / ls() ls(pat = "e") ls(pat = "^m") ls.str() rm(list = ls()) • #列出 当前工作空间中的所有对象 #列出 名字包含e的所有对象 #列出 名字以m开头的所有对象 #列出 当前所有对象的详细信息 #删除 (几乎)一切对象 注意:由于操作函数等是依附于数据的，所以以上操作给出的结果是数据，或是自定义的 函数。 输入 输入 • 任何输入都要考虑如下几个问题: • 从哪儿输入? • 手工 • 文件 • 数据库 • 怎么输入? • 对应的语句是什么?选项是什么? • 输入到哪里去? 输入(手工输入) • 直接在R里面手工输入 x <- 1 x1 <- 2 -> x2
y = c(1, 2, 3, 4)
z = matrix(1:6, 2, 3)
• 变量的命名规则: • 以字符开头
• 可以由字母、下划线“_”、点“.”、数字等组成
• 大小写敏感(例:A与a不能通用)
• 常量名pi，特殊字符NA，NaN 等已被系统保留
• 为避免混淆，尽量避免使用系统已经保留的变量名，函数名，等等。

输入(文件输入)
• 文件就是在计算机上数据的载体，数据存放在文件中 • 两类文件:
• ASCII文件(American Standard Code for Information Interchange， 美国标准信息交换代码)，有时候就简称为文本文件。
• 常见的有.txt文件，.csv文件，等等。
• 优点是在任何机器上都可以读出来，缺点是比较占空间。
• 二进制文件(所有非ASCII文件)
• Excel文件，jpg文件，eps文件，……
• 优点:省空间，除了内容之外还有诸如色彩，排版等其他信息
• 缺点:必须知道内部数据如何存储的才能正确读取
• 注意:一个文件到底是ASCII文件还是二进制文件，不能单凭后缀名来判断。

输入(文本文件输入)
• 问题1:文件放在哪里?(路径) • 问题2:文件名称是什么?
• 问题3:如何读取?
# set working directory
work_dir = “D:/FDUCloud/StatSoft/03” setwd(work_dir)
# input veg data set, method 1
data1 = read.table(file = “veg.txt”, header = T) data1
# input veg data set, method 2
data2 = read.csv(file = “veg.csv”, header = T) data2

输入(文本文件输入)
• 其实这些方法功能比你们想象的强大得多:
read.table(file, header = FALSE, sep = “”, quote = “\”‘”,
dec = “.”, numerals = c(“allow.loss”, “warn.loss”, “no.loss”), row.names, col.names, as.is = !stringsAsFactors,
na.strings = “NA”, colClasses = NA, nrows = -1,
skip = 0, check.names = TRUE, fill = !blank.lines.skip, strip.white = FALSE, blank.lines.skip = TRUE,
comment.char = “#”,
allowEscapes = FALSE, flush = FALSE,
stringsAsFactors = default.stringsAsFactors(),
fileEncoding = “”, encoding = “unknown”, text, skipNul = FALSE)

输入(文本文件输入)
•
•
• •
•
•
read.table
(file, header = FALSE, sep = “”, …)
file
:文件名(包含路径),字符串形式，这个文件名指向你电脑上某个文件的路径
header
:用来表示文件第一行是否含有表头(列名)
sep:
分隔符(separator)，它以字符串形式表示各列是如何分割的。比方说，如
果你有一个用逗号分隔的文件，那么分隔符就是一个逗号。有时文件是用分号、制
表符或是空格分隔的，你就可能需要告诉
read.table()函数分隔符是什么。
Colclasses
: 是一个字符向量，它的长度等于数据集的列数，指明了数据集里每
一列的数据类型。例如，第一列是数值型数据，第二列是逻辑性数据，第三列是因
子。ColClasses
并不是一个必需的向量，但它能告诉read.table()函数每一列的
数据类型。
nrows
:表示数据集的行数，它也不是必需的，但你可以用它来指定行数

输入(文本文件输入)
• read.table(file, header = FALSE, sep = “”, …) •
• •
comment.char:是一个字符串，表示注释标记。一般默认是#(井号)。在注释标
记右边的所有内容都会被忽略。你可以指定其它符号作为注释符号，那么文件中以
注释符号开头的行都会被忽略掉。
skip:表示从第一行开始跳过一定的行数。因为有时在文件的开头会有一些表头信
息或是一些非数据的区域，如果你想直接跳过你可以告诉read.table()函数来跳
过这一部分，例如前十行或是前一百行，只读取那之后的数据。
stringAsFactor: 它的值默认为“真”，通过它可以选择是否把字符变量编码为
因子。默认情况下，只要read.table()函数读到一列字符变量，那么它就会认为
这个变量是一个因子变量。但是如果你不想把这个变量当作因子变量的话，那么你
就可以把stringAsFactor设置为假。

输入(文本文件输入)
• 其实这些方法功能比你们想象的强大得多:
read.csv
(file, header = TRUE, sep = “,”, quote = “\””,
dec = “.”, fill = TRUE, comment.char = “”, …)
read.csv2
(file, header = TRUE, sep = “;”, quote = “\””,
dec = “,”, fill = TRUE, comment.char = “”, …)
read.delim
(file, header = TRUE, sep = “\t”, quote = “\””,
dec = “.”, fill = TRUE, comment.char = “”, …)
read.delim2
(file, header = TRUE, sep = “\t”, quote = “\””,
dec = “,”, fill = TRUE, comment.char = “”, …)

输入(二进制文件输入)
• 问题4:什么格式? • Excel
• XML
• 网页
• SPSS，SAS，Stata，……
library(“readxl”)
work_dir = “D:/FDUCloud/StatSoft/03” setwd(work_dir)
data3 = read_excel(path = “veg.xlsx”, sheet = “veg”)

输入(二进制文件输入)
• 其实这些方法功能比你们想象的强大得多:
read_excel
(path, sheet = NULL, range = NULL, col_names = TRUE,
col_types
= NULL, na = “”, trim_ws = TRUE, skip = 0, n_max = Inf,
guess_max
= min(1000, n_max))
read_xls
(path, sheet = NULL, range = NULL, col_names = TRUE,
col_types
= NULL, na = “”, trim_ws = TRUE, skip = 0, n_max = Inf,
guess_max
= min(1000, n_max))
read_xlsx
(path, sheet = NULL, range = NULL, col_names = TRUE,
col_types
= NULL, na = “”, trim_ws = TRUE, skip = 0, n_max = Inf,
guess_max
= min(1000, n_max))

输入(数据库输入)
• 实际商业操作中，数据都是放在数据库中的
• 常见的数据库有SQLServer,MySQL,PostgreSQL，DB2，SyBase，
Oracle，等等。
• 每种数据库都有自己的实现方式
• 为了便于后台的数据库与前台的应用软件之间的数据交换，也便于不 同数据库之间相互交换数据，规定了“接口”，常见的有两种:ODBC， JDBC。
• R中提供了RODBC和RJDBC两种包
• 其他的包，还有RMySQL，ROracle，RPostgreSQL，RSQLite等等。
• 我们不讲这些包具体怎么用，但是后面会讲SQL语言。

输出

输出
• 任何输出都要考虑如下几个问题:
• 输出什么? • 数据
• 基于数据的统计(后面“描述统计”讲) • 数据分析的结果(后面各个模型章节讲)
• 以什么形式输出? • 数据本身
• 图形
• 输出到哪? • 屏幕
• 文件
• 数据库(不讲)

输出(输出数据)
• 输出到屏幕
(:,,)
( ( ),2,3)
()
(
matrix
matrix
y=
1
6
2
3
z=
rnorm
6
y
print
y
print
round
(z, 2))
• 输出到文件
write.table(z, file = “z.txt”) write.csv(z, file = “z.csv”) save(y, z, file = “yz.RData”) library(“writexl”)
write_xlsx(data.frame(z), path = “z2.xlsx”)
备注:save对应的是load(file = “yz.RData”)

输出(输出数据)
• 与输入类似，输出函数也提供了丰富的细节功能 write.table(x, file = “”, append = FALSE, quote = TRUE, sep = ” “,
eol = “\n”, na = “NA”, dec = “.”, row.names = TRUE, col.names = TRUE, qmethod = c(“escape”, “double”), fileEncoding = “”)
write.csv(…) write.csv2(…)
write_xlsx(x, path = tempfile(fileext = “.xlsx”), col_names = TRUE, format_headers = TRUE)

输出(输出图形)
• 输出到屏幕 x = 1:10
y = x + rnorm(10)
plot(x, y, type = “b”, main = “An Example”)
• 输出到文件
• 文件→另存为→选择对应的图形 • 命令法:
png(“filename.png”)
plot(x, y, type = “b”, main = “An Example”) dev.off()

输出(输出图形)
• R可以输出多种格式的图形，常见的有如下这些: pdf(“filename.pdf”)
win.metafile(“filename.wmf”) png(“filename.png”) jpeg(“filename.jpg”) bmp(“filename.bmp”) postscript(“filename.eps”)
• 点阵图形与矢量图形

输出(输出图形)
• 这些输出图形也有丰富的选项可供定制细节:
pdf(file
= ifelse(onefile, “Rplots.pdf”, “Rplot%03d.pdf”),
width, height,
onefile, family, title, fonts, version,
paper, encoding,
bg, fg, pointsize, pagecentre, colormodel,
useDingbats
, useKerning, fillOddEven, compress)
bmp
filename
“Rplot%03d.bmp”
width
=
480, height
=
480, units
px
pointsize
bg
=
“white”
, res
=
NA
,
family
=
“”
restoreConsole
type
“windows”,
cairo
antialias
(=,
=””, =,
,=, =c( ” “), )
(=,
= = “px”, pointsize = 12,
,
=, =, =TRUE,
=c( ,” “), )
12
TRUE
jpeg
filename
“Rplot%03d.jpg”
width
= =
480, height
480, units
quality
75
bg
=
“white”
, res
NA
family
“”
restoreConsole
type
“windows”
cairo
antialias

输出多种结果到多个地方
sink(“code_2_sub.txt”, append = TRUE, split = TRUE) png(“xy2.png”)
source(“code_2_sub.R”)
sink()
dev.off()
• code_2_sub.R
x = 1:10
y = x + rnorm(10)
print(x)
print(y)
plot(x, y)
• sink函数指定输出的文本文件，append选项为在现有文件后面添加，split 选项表示是否同时有多个输出(即同时在屏幕和文件内输出)。
• sink函数仅对文本结果有效，对于图形结果无效。 • 做完之后别忘了sink(),dev.off()

数据操作 R中常见的数据类型与操作

R语言里面有哪些数据类型?
• 常见的单个数据(类型):
• 逻辑型(布尔型，TRUE/FALSE)
• (实)数值型
• 整数型和双精度浮点数型两种
• 日期型与时间型两种特殊的数值型
• 复数型(本课程不讨论)
• 字符型(单个字符和字符串都属于字符型) • 原生型(本课程不讨论)

数据操作 R中常见的数据类型与操作 逻辑型数据

逻辑型数据的相关操作符
• 逻辑运算(与、或、非)
x1 = c(TRUE, TRUE, FALSE, FALSE)
x2 = c(TRUE, FALSE, TRUE, FALSE)
x1 & x2
x1 | x2
!x1
x1 == x2
xor(x1, x2)

逻辑型数据的相关操作符
• 逻辑运算(存在与所有，TRUE的位置)
x1 = c(TRUE, TRUE, FALSE)
any(x1)
all(x1)
which(x1)
x2 = c(TRUE, TRUE, TRUE)
any(x2)
all(x2)
which(x2)
x3 = c(FALSE, FALSE, FALSE)
any(x3)
all(x3)
which(x3)

逻辑型数据的相关操作符
构造逻辑型数据
x = c(-2, -1, 0, 1, 2) x == 0
x != 0
x<0 x <= 0 x>0 x >= 0

数据操作 R中常见的数据类型与操作 数值型数据

数值型数据处理操作符

数值型数据处理函数
• 既然要处理数据，当然必须有一些数学函数

数值型数据处理函数

数据操作 R中常见的数据类型与操作 补充:概率统计运算

概率统计运算
• 统计软件当然必须有一些统计函数

概率统计运算
• 生成服从标准正态分布的随机数，计算pdf,cdf及其逆函数。 # 生成100个标准正态的随机数
x1 = rnorm(100)
# 生成100个均值为2，标准差为0.5的标准正态随机数 x2 = rnorm(100, mean = 2, sd = 0.5)
x2 = rnorm(100) * 0.5 + 2
# 计算标准正态分布在0.5处的密度函数
dnorm(0.5)
# 计算标准正态分布在0,1.65,1.96处的累计密度函数 pnorm(c(0, 1.65, 1.96))
# 计算标准正态分布0.5, 0.90, 0.95, 0.975分位数 qnorm(c(0.5, 0.90, 0.95, 0.975))
• 注:其他分布请参看R的帮助以及其他的文档。

概率统计运算
• 画标准正态分布密度函数图
x=
seq(from = -3, to = 3, by = 0.1)
y1 =
dnorm(x)
y2 =
pnorm(x)
plot
(x, y1, type = “l”, col = “blue”)
lines
(x, y2, type = “l”, col = “red”)

概率统计运算
• 随机数种子 rnorm(5)
rnorm(5)
set.seed(11) rnorm(5) set.seed(11) rnorm(5)

概率统计运算(for 普通青年) • 有哪些内置的概率分布

概率统计运算(for 文艺青年)
• 如果某些分布R里面没有提供，可以通过他们之间的关系来间接生成

概率统计运算(for 其他青年)

概率统计运算
• 生成多元密度随机数
• 比如生成100个二元正态分布随机数，均值为 1 ，方差-协方差矩阵为
1 0.5 2 0.5 2
• 注意:需要调用MASS包。 library(“MASS”)
mu = matrix(c(1, 2), 2, 1)
sigma = matrix(c(1, 0.5, 0.5, 2), 2, 2) x = mvrnorm(100, mu, sigma)

数据操作 R中常见的数据类型与操作 日期型数据

日期型数据处理
• 日期值通常以字符串的形式输入到R中，然后转化为以数值形式存储的日期变量。 • 函数as.Date()用于执行这种转化。
•
其语法为as.Date(x, “input_format”)，其中x是字符型数据，input_format
则给出了用于读入日期的适当格式。

日期型数据处理
x1 = “2/3/2020”
x1
class(x1)
x2 = as.Date(x1, “%m/%d/%Y”)
x2
class(x2)
# print system date
Sys.Date()
date()
# Woring with formats
today <- Sys.Date() format(today, format="%B %d %Y") format(today, format="%A") 日期型数据处理 • • R的内部在存储日期时，是使用自1970年1月1日以来的天数表示的，更早的日期则 表示为负数。 可以在日期值上执行算术运算 x< - as.Date "1970 01 01" as.numeric as.numeric as.numeric (--) (x) ((-- ((-- as.Date as.Date "1970 "1969 01 12 02")) 02")) # Calculations with with dates startdate <- as.Date("2004-02-13") enddate <- as.Date("2009-06-22") enddate - startdate # Date functions and formatted printing today <- Sys.Date() dob <- as.Date("1956-10-12") difftime(today, dob, units="weeks") 日期型数据处理 • R语言中用了特殊的数据类型来表示日期和时间 • • • • • • • • 日期用 Date类来表示 日期不包括时间，只表示某年某月某日，可以用年月日这样的格式来表示。 内部日期是以 1970年1月1日至今的天数来存储 内部时间则是以 1970年1月1日至今的秒数来存储。 日期作为 Date类对象来存储，在R里的工作方式是输入一个字符串, 比如1970- 01 -01，然后用as.Date()函数转换成日期。 时间则是由两个不同的类: POSIXct POSIXlt 数据操作 R中常见的数据类型与操作 时间型数据 时间型数据处理 将字符串转换为时间型数据 有两个字符串:“October 7, 2015 10:40”, “December 9, 2011 9:10 ”,如果把这些字符串转 换成时间对象，即可以使用strptime()函数，其中%B代表月份，%d代表日，%Y是四位数字的年份， %H是小时，%M是分钟。 • datestring = c("October 7, 2015 10:40", "December 9, 2011 9:10") datestring class(datestring) x = strptime(datestring, "%B %d, %Y %H:%M") x Sys.getlocale() Sys.setlocale("LC_ALL", "English") x = strptime(datestring, "%B %d, %Y %H:%M") x class(x) Sys.setlocale("LC_ALL", "Chinese") 时间型数据处理 • 时间则是用两种类型来表示，一种是POSIXct，另一种是POSIXlt。 • • • • POSIX是一个计算标准家族，它规定了特定类型的计算机处理事件和数 据表示的方法。 在POSIXct类里，时间是用非常大的整数来表示，如果想要把时间存在 数据框里的话，这种类非常有用，因为它基本上就是一个很大的整数向 量。 而POSIXlt实际上把时间当作列表来存储，它还存在一些跟给定时间相 关的信息，例如这个时间是星期几，是一年中的第几天，是几号?是几 月? 使用as.POSIXct()或as.POSIXlt()函数，可以让对象在POSIXct和 POSIXlt之间来回类型转换 时间型数据处理 两种时间型数据的结构 # POSIXct时间型数据的结构 x1 <- Sys.time() x1 class(x1) str(x1) unclass(x1) str(unclass(x1)) names(x1) # 转换成POSIXlt时间型数据及其结构 x2 <- as.POSIXlt(x1) x2 class(x2) str(x2) unclass(x2) str(unclass(x2)) names(x2) 时间型数据处理 • 对于POSIXct格式，仍可以用Sys.time()函数来得到当前时间; • • • 如果读POSIXct对象去类型化，得到一个非常大的整数，即自1970年1月1日 至今的秒数; 如果对这个对象使用列表的$操作符的话，会发现得到的是一个错误，因为 POSIXct类的对象没有这样的列表元素; 如果想得到这些列表元素，必须用as.POSIXlt()函数把它转换成POSIXlt类， 才可以得到秒数。 x1 = Sys.time() x2 = as.POSIXlt(x1) x1$sec # 本句执行时会报错 x2$sec 时间型数据处理 • • 注意:不能把不同的类混淆在一起。 有了日期或时间格式的数据，可以很方便对它们进行运算，比方说可以 对日期做减法，比较日期大小等。 Sys.setlocale ("LC_ALL", "English") x = as.Date ("2015-10-07") y = strptime ("9 Jan 2011 11:34:21", "%d %b %Y %H:%M:%S") class( x) class( y) x-y # 本句执行时会报错 x = as.POSIXlt (x) x x-y Sys.setlocale ("LC_ALL", "Chinese") 时间型数据处理 • • R里日期和时间运算的好处在于它们能自动处理好一些问题，例如闰年、闰秒、夏 令时以及时区等。 以下以2020年为例，2020是闰年，所以它有2月29号，3月1日和2月28之间的差实 际上是2天，而不是像平年那样相差1天。 x <- as.Date("2020-03-01") y <- as.Date("2020-02-28") x–y x <- as.Date("2019-03-01") y <- as.Date("2019-02-28") x-y 时间型数据处理 • 同样，可以取北京时间x,以及格林尼治时间(Greenwich Mean Time)y,然后求差 值。虽然看起来一样，实际上相差了8小时，因为这里有时区的改变。(我们在东8 区) x <- as.POSIXct("2020-03-10 01:00:00",tz="PRC") y <- as.POSIXct("2020-03-10 01:00:00",tz="GMT") x-y x <- as.POSIXlt("2020-03-10 01:00:00",tz="PRC") y <- as.POSIXlt("2020-03-10 01:00:00",tz="GMT") x-y x <- as.POSIXct("2020-03-10 01:00:00",tz="PRC") y <- as.POSIXlt("2020-03-10 01:00:00",tz="GMT") x-y 数据操作 R中常见的数据类型与操作 字符型数据 字符型数据处理函数 字符型数据处理函数 字符型数据处理函数 字符型数据处理函数 • 字符型转成因子型(分成有序和无序两种)。 • 什么时候需要转成因子型? x = c("Type C", "Type B", "Type B", "Type A", "Type A", "Type A") class(x) str(x) x1 = factor(x) x2 = factor(x, ordered = T) x3 = factor(x, ordered = T, levels = c("Type C", "Type B", "Type A")) class(x1); class(x2); class(x3) str(x1); str(x2); str(x3) 数据操作 R中常见的数据类型与操作 集合运算 集合运算 • 对各种类型数据都适用，此处仅以数值型数据为例。 x1 = c(1, 2, 3) x2 = c(3, 4, 5) union(x1, x2) intersect(x1, x2) setdiff(x1, x2) setequal(x1, x2) x1 %in% c(1, 2, 3, 4) x1 %in% x2 数据操作 R中常见的数据结构与操作 R语言里面有哪些数据结构? • 常见的单个数(据): • 标量 • 从单个类型的数(据)拓展开去: • 向量(一维数组) • 矩阵(二维数组) • 高维数组 • 不同类型的组合: • 数据框(data.frame) • 列表 • 注意对比CS中数据结构 • 所有这些在R里都叫做“对象”，R中的对象还包括函数，命令，等等。 “一切皆对象” 数据操作 R中常见的数据结构与操作 标量 标量 • 标量即只有一个数(数值，日期，时间，字符串，逻辑，等等) • 标量可以看成是一个长度为1的向量(数组) a1 = 3 a1 class(a1) a2 <- "Hello, World!" a2 class(a2) a3 <- TRUE -> a4
a3
a4
class(a3)
class(a4)
a5 = as.Date(“2019-10-15”)
a5
class(a5)

数据操作 R中常见的数据结构与操作 向量

向量
• 向量就是所有元素类型相同的一维数组 • 常见向量(以及向量赋值运算)
a1 = c(1, 2, 5, 3, 6, -2, 4)
a2 = c(“one”, “two”, “three”)
a3 = c(TRUE, TRUE, TRUE, FALSE, TRUE, FALSE)
a4 = 1:10
a5 = seq(from = -5, to = 5, by = 3)
a6 = rep(2, 5)
•
访问向量中的元素
a1[3]
a1[c(1,3,5)]
a1[2:6]
• c函数:Combine Values into a Vector or List

向量
• 注意:向量中所有的元素都必须是同一个类型 • 我们试一下以下操作
d1 = c(1, “one”)
d2 = c(1, TRUE)
d3 = c(“one”, TRUE)
d4 = c(1, “one”, TRUE)
• 我们在数据整理部分会专门再讲数据类型转换

数据操作 R中常见的数据结构与操作 矩阵

矩阵
• 创建矩阵(所有元素类型相同的二维数组) y = matrix(1:20, nrow = 5, ncol = 4)
y
cells = c(1, 26, 24, 68)
rnames = c(“R1”, “R2”)
cnames = c(“C1”, “C2”)
mymatrix = matrix(cells, nrow = 2, ncol = 2, byrow = FALSE,
dimnames = list(rnames, cnames))
mymatrix

矩阵
• 矩阵元素提取
x = matrix(1:10, nrow = 2)
# print x
x
# print the 2nd row of x
x[2, ]
# print the 2nd column of x
x[, 2]
# print the element on the 1st row and 4th column of x
x[1, 4]
# print the element on the 1st row, 4th and 5th column of x
x[1, c(4, 5)]
# print the element on the 1st row, 4th column of x (twice)
x[1, c(4, 4)]

矩阵
• 对矩阵的计算
x = matrix(1:12, 3, 4)
x
t(x)
sqrt(x)
mean(x)
sum(x)
colSums(x)
colMeans(x)
rowSums(x)
rowMeans(x)
apply(x, 1, max)
apply(x, 2, min)

数据操作 R中常见的数据结构与操作 多维数组

多维数组
• 创建多维数组(所有元素类型相同) dim1 = c(“A1”, “A2”)
dim2 = c(“B1”, “B2”, “B3”)
dim3 = c(“C1”, “C2”, “C3”, “C4”) z = array(1:24, c(2, 3, 4),
dimnames = list(dim1, dim2, dim3))
z
• 提取元素 z[1, 2, 3]

为什么要做随机模拟

回顾上学期内容
• 𝑋 ∼ 𝜒2(𝑑)，对应的概率密度函数: 1 𝑑−1 −𝑥/2
• 2𝑑Γ(𝑑)𝑥2 𝑒 2
• 𝑋 ∼ 𝑡(𝑑) ，对应的概率密度函数: Γ 𝑑+1 2 −𝑑+1
•21+𝑥2
𝑑𝜋Γ 𝑑
• 𝑋 ∼ 𝐹(𝑑1, 𝑑2) ，对应的概率密度函数:
𝑑1𝑥 𝑑1𝑑𝑑2 2
𝑑1𝑥+𝑑2 𝑑1+𝑑2
𝑥𝐵(𝑑1,𝑑2) 22
• 看起来很烦?其实这些都是简单的!因为至少你还能把数学表达式给写 出来!!
𝑑 2
•

不信?你试试这几个!
•𝑋∼𝑁(𝜇 ,𝜎2)，Y∼𝑁(𝜇 ,𝜎2) 𝑥𝑥 𝑦𝑦
• 请写出𝑋/𝑌，𝑋2𝑌， 𝑋𝑌分布的p.d.f.
• 这些理论推导即便对于概率论专家来说也不是一件简单的事情，而且还
不一定写得出来!
• 实际中碰到这些问题怎么办?有没有简单点的方法，哪怕不是那么精确， 但是至少能够给我们一个答案?

不信?你试试这几个!
•𝑋∼𝑁(𝜇 ,𝜎2)，Y∼𝑁(𝜇 ,𝜎2) 𝑥𝑥 𝑦𝑦
• 请写出𝑋/𝑌，𝑋2𝑌， 𝑋𝑌分布的p.d.f.，再求期望、方差……
• “简单”“粗暴”的模拟法
• 生成1,0000个𝑋，再生成1,0000个Y
• 计算上面的数值，然后画直方图，算期望和方差。 • 1,0000不够就10,0000，不行再加
• 这是没有办法的办法，但总比没有办法强! • 回想一下蒲丰投针问题

为什么要做随机模拟
• 在很多问题中，人们需要计算一些牵涉到很复杂的随机变量的 期望、方差等问题
• 由于这些随机变量形式复杂，直接计算往往比较困难
• 可以通过随机模拟的方法，将这类概率问题“转化”成统计问
题:将感兴趣的随机变量当作一个“总体”，从该总体中抽取 样本，基于这些样本数据来估计需要计算的量。

随机模拟的应用
• 随机模拟是一种需要依靠计算机实现的技术，用于处理与随机 变量相关的问题，其特点是可以模拟较为复杂的不确定性环境，
为管理人员提供决策支持。
• 随机模拟有众多应用: • 金融衍生品的策略/建模 • AlphaGo/Alpha Zero
• 传染病传染过程模拟 •……

随机模拟案例1:Conley水产公司

Conley水产公司
• 例11.1:Conley水产公司的经 营
• 每天能捕获3500磅鳕鱼，成本 为$1,0000;
• 港口1
• 价格 𝑃 = $3.25/磅
• 需求量𝐷总是大于3500磅
• 港口2
• 价格𝑃 ∼ 𝑁(3.65,0.22)
• 需求量𝐷 服从一个离散分布
• 价格与需求量相互独立
需求量可能的取值 (磅) 概率函数
0 0.02
1,000 0.03
2,000 0.05
3,000 0.08
4,000 0.33
5,000 0.29
6,000 0.20

Conley水产公司
• 管理问题:应该去哪一个港口销售?
• 港口1
• 收入 I = 3.25×3500 = 11375
• 利润 F=3.25×3500−10000=1375
• 港口2
• 收入I = ቊP × 3500，若D ≥ 3500 = P × min(3500,D)
P × D，若D < 3500 • 利润F=P×min 3500,D −10000 • 显然去港口2销售的利润𝐹是个随机变量，为了能做决策，需要 了解𝐹的分布特征 • 如何来求𝐹的概率分布? Conley水产公司 • 假设Conley公司财大气粗，不在乎盈亏，那么可以做这样一 个实验:在接下去的200天里，每天去捕鱼，然后去港口2销售， 如此就可以获得总体F的一个容量为200的样本(实现值)，接 下来用这200个样本数据来估计F的分布特征，比如期望、标准 差，等等 • 不过，如果Conley公司经费没有那么充沛，应该没有闲情雅 致去这样捕鱼...... • 事实上，这个思想未必要通过真实的实验来达成，Conley公 司可以借助随机模拟来实现。 Conley水产公司 • 第1步:按照𝐷的概率分布，生成𝑚个(比如𝑚 = 10000)服从 该分布的随机数(样本实现值) • 𝑑1,⋯,𝑑𝑚 • 第2步:按照𝑃的概率分布(即𝑁(3.65,0.22))，生成𝑚个服从 该分布的随机数 •𝑝,⋯,𝑝 1𝑚 • 第3步:根据计算公式，计算𝑚个利润的实现值 • 𝑓 =𝑝 min 3500,𝑑 −10000，𝑖=1,⋯,𝑚 𝑖𝑖𝑖 • 第4步:将算得的𝑚个利润实现值当作𝑚个样本观测值，应用统 计方法加以分析并回答问题。 Conley水产公司 # set the #(simulation) m = 50000 # generate demand probability = c(0.02, 0.03, 0.05, 0.08, 0.33, 0.29, 0.20) income = seq(0, 6000, by = 1000) d = sample(income, m, replace = T, prob = probability) d[d > 3500] = 3500
# generate price
p = rnorm(m, mean = 3.65, sd = 0.2)
# calculate revenue, f
f = p * d – 10000
# calculate summary statistics of f
hist(f); mean(f); sd(f); quantile(f,0.25)

随机模拟案例2:小饭店

小饭店
• 假定你要开一个小饭店 • 月度固定成本
• 租金等固定成本:𝑈 = 3995
• 月度可变成本 • 劳动力成本:𝐿
• 月度膳食销售数量:𝑀
• 月度膳食可变成本:𝐹 = 11 × 𝑀
• 月度总成本
• 𝐿+𝑈+𝐹=𝐿+𝑈+11×𝑀

小饭店
• 假定你要开一个小饭店
• 月度收入
• 每餐饭的价格:𝑃
• 月度收入:𝑅=𝑃×𝑀
• 因此，月度利润为
• 𝑋=𝑃×𝑀− 𝐿+𝑈+11×𝑀 = 𝑃−11 ×𝑀−𝐿−3995
• 哪些量是随机变量?
• 劳动力成本:𝐿
• 月度销售膳食的数量: 𝑀
• 每餐饭的价格: 𝑃
• 月度利润:𝑋( 𝑋为随机变量的函数，也是随机变量)

小饭店
• 劳动力成本:L
• 假设劳动力成本会分布在每月5040至6860之间，在这个范围内服从一个均匀分
布
• 月度膳食销售数量: M
• 假设每月销售的膳食数量服从一个均值为3000、标准差为1000的正态分布 • 每餐饭的价格:P
• 假设每餐饭的价格服从一离散型概率分布 • 假定L，M，P独立。
场景 每餐饭的价格 概率
非常兴旺的市场 $20.00 0.25
兴旺的市场 $18.50 0.35
不是非常兴旺的市场 $16.50 0.30
不兴旺的市场 $15.00 0.10

小饭店
• 感兴趣的随机变量是:月度利润𝑋
• 计算出𝑋的概率分布比较麻烦。然而，𝑋是𝐿、 𝑀与𝑃的
函数，其随机性完全由这几个数决定。
•要准确地理解地分布，对𝐿、 𝑀与𝑃 所做的假设很关键。
• 在实际应用中，需要考虑
• 模型与实际情况的吻合度如何?
• 假设与实际情况有多大偏差? • 分布的假设
• 独立性的假设

小饭店
• 感兴趣的问题
• 𝑋服从什么概率分布?
• 𝑋的期望与标准差分别为多少? • 𝑃(𝑋 ≤ 5000)
• 𝑃(𝑋 ≥ 6667)
• 其他感兴趣的问题?

小饭店
生成n组随机数
(如n = 2000)
(L , M , P ) 111
(L , M , P ) 222

(L , M , P ) nnn
X = (P − 11)  M − L − 3995
– E(X)
– SD(X)
– P(X  $5,000) – P(X  $6,667) – 其他量 的估计值
X =(P−11)M −L−3995 1111
X =(P −11)M −L −3995 2222

X =(P −11)M −L −3995 nnnn
用生成的随 机数来估计 感兴趣的量

小饭店
◼ 月底利润计算的2000次模拟数据:
模拟
次数
月度 月度劳动 利润 力成本
月度膳食 每餐饭 销售数量 的价格
1 2 3
-2,067.54 5,241.67
20,073.45 13,484.71 10,843.88
6,592.06 6,781.13 5,459.35
4,088 18.50 3,235 18.50 3,691 16.50
956 18.50
13,278.86 5,263.36
2,504 20.00
4 5 6 7 …
1,597.77 5,200.25
1,962 16.50
-3,592.04 6,561.69
1,741 15.00
………… 15,106.93 6,566.22 3,422 18.50
◼ 基于这些模拟的月度利润额，我们希望: ◼ 估计 X 的概率分布
◼ 估计 X 的期望与标准差
◼ 估计 P(X  $5000)
◼ 估计 P(X  $6667)
……
……
… 1999 2000
12,777.25 6,500.11
4,231 16.50

小饭店
•基于𝑋的模拟数据画直方图，可以作为𝑋 概率分布的直 观估计

小饭店
◼基于2000个模拟数据 (当作样本) ◼ 估计期望𝐸 𝑋 = 10854.66
◼ 估计标准差𝑆𝐷 𝑋 ◼ 估计𝑃 𝑋 ≤ 5000 ◼ 估计𝑃 𝑋 ≥ 6667
= 8802.35 = 0.2845 = 0.6405
◼ 代码:自己试着做!

小饭店
• 除了考虑月度利润，或许分析饭店的年度利润也有价值
• 完全可以用随机模拟的方法来分析年度利润，但需要给出更多
的假设
• 保留月度利润分析时所做的所有假设，并进一步假设:
• 对于一年中的每个月均引入一个随机变量描述当月膳食销售数量，记为 𝑀 , ⋯ , 𝑀 ，假设这些随机变量独立同分布
1 12
• 假设月度每餐饭的价格𝑃在一整年中保持不变
• 假设月度劳动力成本𝐿在一整年中保持不变

小饭店
• 根据𝑀 , ⋯ , 𝑀 ，𝑃与𝐿各自的分布独立地产生随机数，由此模 1 12
拟出每个月的利润额
• 𝑋 = 𝑃 − 11 × 𝑀 − 𝐿 − 3995 𝑖𝑖
• 由此年度利润额为 •𝐴=𝑋+⋯+𝑋
1 12
• 仍然需要估计: • 𝐴的概率分布
• 𝐴的期望与标准差 • 概率𝑃(𝐴 ≤ 60000) • 概率𝑃(𝐴 ≥ 60000)

小饭店

随机模拟案例3 大数定律与中心极限定理

大数定律与中心极限定理
• 常见的统计量很多都是均值:
• 已有两组随机数:𝑋 ,𝑋 ,⋯,𝑋 ，𝑌 ,𝑌 ,⋯,𝑌
• 均值: • 方差:
ത 𝑋=
𝑋 +⋯+𝑋 1𝑛
𝑛
12𝑛12𝑛
ത2 ത2 𝑋−𝑋 +⋯+𝑋−𝑋
𝐶𝑜𝑣 𝑋,𝑌 =
11 𝑛𝑛
𝑉𝑎𝑟 𝑋 =
1𝑛
𝑛
• 协方差:
• 相关系数等:均可以看做是均值的函数
തത തത (𝑋 −𝑋)(𝑌 −𝑌)+⋯+(𝑋 −𝑋)(𝑌 −𝑌)
𝑛

大数定律与中心极限定理
均值很忙!
均值很重要! 均值好不好?
均值好在哪?

大数定律与中心极限定理
• 以下，我们讨论范围仅限于独立同分布数据。 • 假定 𝑋 , ⋯ , 𝑋 为一组独立同分布的随机变量
ത •𝑋𝑛=
• 考虑:
𝑋 +𝑋 +⋯+𝑋 12𝑛
𝑛
1𝑛
𝑋 ,⋯,𝑋 1𝑛
独立 随机
ത 𝑋是个随机变量
തത
•问:𝑋 的均值是什么?𝑋 分布是什么?

随机模拟案例3.1 大数定律

大数定律
 Law of Large Number (L.L.N.) 𝑋 ,𝑋 ,⋯,𝑋 独立同分布
12𝑛
 𝑋𝑖 所服从的分布均值为 𝜇
𝑃
 随着𝑛 → ∞，有𝑋(𝑛) → 𝜇
 随着样本量的增大，样本均值趋向于总体均值
 解释:
 抛硬币，看正反面(扔多了，正反一半对一半)  赌博，玩多了总是输的次数多
 很长的一个时间段内，股价涨跌一半对一半。
ത

大数定律
 生成10000个独立同分布的随机 变量𝑋 ,𝑋 ,⋯,𝑋 ，服从参
1 2 10000
数为0.5的Bernoulli分布  计算前n个数的算术均值
ത 𝑋 +𝑋 +⋯+𝑋
𝑋(𝑛)=12 𝑛 𝑛
 重复4次，观察现象

大数定律
 生成10000个独立同分布的随机 变量𝑋 ,𝑋 ,⋯,𝑋 ，服从参
1 2 10000
数为0.5的Bernoulli分布  计算前n个数的算术均值
ത 𝑋 +𝑋 +⋯+𝑋
𝑋(𝑛)=12 𝑛 𝑛
 重复500次，观察

大数定律
 生成10000个独立同分布的随机 变量𝑋 ,𝑋 ,⋯,𝑋 ，服从
1 2 10000 [0,1]上均匀分布
 计算前n个数的算术均值 ത 𝑋 +𝑋 +⋯+𝑋
𝑋(𝑛)=12 𝑛 𝑛

大数定律
 生成10000个独立同分布的随机 变量𝑋 ,𝑋 ,⋯,𝑋 ，服从标
1 2 10000
准正态分布
 计算前n个数的算术均值
ത 𝑋 +𝑋 +⋯+𝑋
𝑋(𝑛)=12 𝑛 𝑛

随机模拟案例3.2 中心极限定理

中心极限定理
 生成10000个独立同分布的随机 变量𝑋 ,𝑋 ,⋯,𝑋 ，服从参
1 2 10000
数为0.5的Bernoulli分布
 我们画出𝑛为100,400,10000时
ത ，𝑋(𝑛)的直方图

中心极限定理
 生成10000个独立同分布的随机 变量𝑋 ,𝑋 ,⋯,𝑋 ，服从
1 2 10000 [0,1]上均匀分布
 我们画出𝑛为100,400,10000时
ത ，𝑋(𝑛)的直方图

中心极限定理
 生成10000个独立同分布的随机 变量𝑋 ,𝑋 ,⋯,𝑋 ，服从标
1 2 10000
准正态分布
 我们画出𝑛为100,400,10000时
ത ，𝑋(𝑛)的直方图

中心极限定理
• 均值不变
• 方差越来越小
ത
• 𝑋(𝑛)的分布对称，且中间高、两头低

中心极限定理
Bernoulli(0.5)
Unif(0,1)
N(0,1)
𝑛
ത
|𝑋 𝑛 −𝜇|
𝑠(𝑛)
ത
|𝑋 𝑛 −𝜇|
𝑠(𝑛)
ത
|𝑋 𝑛 −𝜇|
𝑠(𝑛)
100
0.00098
0.049596
0.001929
0.028187
0.001377
0.097053
400
0.00054
0.025320
0.001214
0.014755
0.000501
0.047908
10000
0.00010
0.005083
0.000050
0.002920
0.000140
0.010287
三个现象:
തത
1. 𝑋 𝑛 −𝜇 →0，即𝑋(𝑛)→𝜇
2. 𝑠(𝑛)→0
3. 𝑠(𝑘2𝑛) ≈ 𝑠(𝑛)/𝑘 (上表中𝑘 = 2，10)

中心极限定理
ത
• 𝑋(𝑛)的分布到底趋于哪个分布?
• 我们引入一个新的工具:QQ图(Q = Quantile)
ത
ത
ത 2 𝑋(𝑛)−𝜇
•如果𝑋𝑛∼𝑁𝜇,𝜎 ，则 𝑛𝑛𝜎
𝑛∼𝑁(0,1)
• 𝑋(𝑛)−𝜇𝑛的𝛼分位数 =标准正态分布的𝛼分位数(Φ−1(𝛼)) 𝜎
𝑛
𝑛
ത −1
•𝑋(𝑛)的𝛼分位数 = 𝜇 +𝜎 Φ (𝛼) 𝑛𝑛
ത
• 如果我们将𝑋(𝑛)和标准正态分布对应的分位数分别作为横坐标与纵坐标，
一一对应画下来，应该正好是一条直线

中心极限定理
 生成10000个独立同分布的随机 变量𝑋 ,𝑋 ,⋯,𝑋 ，服从参
1 2 10000
数为0.5的Bernoulli分布
 我们画出𝑛为100,400,10000时
ത ，𝑋(𝑛)的QQ图

中心极限定理
 生成10000个独立同分布的随机 变量𝑋 ,𝑋 ,⋯,𝑋 ，服从
1 2 10000 [0,1]上均匀分布
 我们画出𝑛为100,400,10000时
ത ，𝑋(𝑛)的QQ图

中心极限定理
 生成10000个独立同分布的随机 变量𝑋 ,𝑋 ,⋯,𝑋 ，服从标
1 2 10000
准正态分布
 我们画出𝑛为100,400,10000时
ത ，𝑋(𝑛)的QQ图

中心极限定理
• 没有假定𝑋 具体是哪个分布 𝑖
• 也没有假定𝑋 是离散还是连续 𝑖
ത
• 算术均值𝑋的极限分布都是正态分布

中心极限定理
• Central Limit Theorem (C.L.T.) • 𝑋 ,𝑋 ,⋯,𝑋 独立同分布
12𝑛
• 𝑋 均值为 𝜇，方差为𝜎2
ത𝑑
𝜎2 • 𝑋 → 𝑁(𝜇, 𝑛 )
• 备注:
• 样本量要足够大，一般要求不小于30
• 常见错误1:只要样本量足够大，样本就趋于正态分布
• 常见错误2:考试成绩被正态!90分以上控制在30%! 

中心极限定理
• 均值好不好?
• 均值好在哪?
• 大数定律:随着样本量的增加，样本均值趋于分布均值
• 中心极限定理:随着样本量的增加，样本均值的分布趋于正态分布(均值与原 分布的均值相同，标准差为原分布的标准差除以开根号样本量，或者说方差为 原分布的方差除以样本量)
• 注意:要求样本独立同分布

随机模拟:小结

小结
• 随机模拟可以完成一些普通理论推导很难完成的不确定性分析 • 随机模拟可以处理复杂的和动态的系统
• 由于随机模拟技术本身具有内在的随机性，因此其给出的结果 并非精确结果 (可以通过增加模拟次数来减少随机性)，在解
释模拟结果时需要注意至一点
• 随机模拟的优点:可以处理任意分布;实现简便
• 随机模拟的限制:不提供最优决策;有效性依赖于分布假设

小结
• 在实际应用中，当需要通过随机模拟技术来帮助进行管理决策 时，需要注意模拟机结果的可靠性至少依赖于以下3点:
• 对于所面对管理问题的背景有深刻的了解
• 对于随机模拟中牵涉到的概率论与统计学基本原理有正确的理解 • 能将管理问题正确地转化成合理的随机模拟模型

模型回归系数评估

模型回归系数评估:回顾一下估计
𝑦=𝛽+𝛽𝑥 01
መመ 𝑦=𝛽+𝛽𝑥
01

模型回归系数评估:回顾一下估计
•𝑦=𝛽+𝛽𝑥+𝜀 𝑖01𝑖𝑖
መመ𝑛መመ2 •最小二乘法:找𝛽,𝛽，使得σ (𝑦−𝛽 −𝛽𝑥)最小
01 𝑖=1𝑖01𝑖
• 最小二乘估计(Least Square Estimator，LSE) መ σ 𝑛 ( 𝑥 𝑖 − 𝑥 ҧ ) ( 𝑦 𝑖 − 𝑦ത )
• 𝛽 = 𝑖=1
1 σ𝑛 𝑥𝑖−𝑥ҧ 2
መመ
• 𝛽 = 𝑦ത − 𝛽 𝑥 ҧ
• 附带估计
01
መመ
• 𝑦ො = 𝛽 + 𝛽 𝑥
𝑖=1
(𝑦 的拟合值) 𝑖01𝑛𝑖2 𝑖
෢2 2 σ𝑖=1𝑦𝑖−𝑦ො𝑖 2
•𝜎 =𝑠 = 𝑛−2 (𝜎的估计)
• 注意:所有的估计都是随机变量!!!

模型回归系数评估
መ
𝛽 点估计
𝑖
𝑖
෢2 𝑠𝑖= 𝜎
假设检验
𝐻:𝛽=0 0𝑖
መ
𝛽 的样本标准差
𝑖
෡
𝛽𝑖−𝛽𝑖服从自由度为𝑛 − 1 − 1的𝑡分布 𝑠𝑖
先看这个结果，我们下面一个个解释。

模型回归系数评估(估计的角度)
• 𝛽 和𝛽 的点估计(见第二列) 01𝑛
መ σ ( 𝑥 𝑖 − 𝑥 ҧ ) ( 𝑦 𝑖 − 𝑦ത )
• 𝛽 = 𝑖=1
1 σ𝑛 𝑥𝑖−𝑥ҧ 2
መመ
• 𝛽 = 𝑦ത − 𝛽 𝑥 ҧ
01
𝑖=1
መመ መመ
• 𝛽 与𝛽 都是随机数(𝛽 与𝛽 未必等于真实的𝛽 与𝛽 )
010101
መመ
• 问1:𝛽 与𝛽 的均值是多少?
01
• 答1:就是𝛽 和𝛽 01
መመ
• 问2:𝛽 与𝛽 的标准(误)差𝜎 与𝜎 是多少?
01 01
• 答2:见第Std. Error列(注意:是估计出来的标准误差𝑠 与𝑠 ) 01

模型回归系数评估(估计的角度)
መመ
• 问3:𝛽 与𝛽 服从什么分布?
01
• 答3:正态分布(前提是𝜀𝑖服从正态分布)
መመ
• 问4:(𝛽 −𝛽 )/𝑠 与(𝛽 −𝛽 )/𝑠 服从什么分布?
000111 • 答4:自由度为𝑛 − 2的𝑡分布

模型回归系数评估(检验的角度)
መመ
• 如果真实的的𝛽 或 𝛽 等于0，但是估计出来的𝛽 与𝛽 由于随机性的影响而不等于0。
•检验3:𝐻:𝛽 =0 01
01መመ 01
我们还是按照估算出来的𝛽 与𝛽 去预测的话，会产生系统性(模型识别)偏差。 01
• 我们要把一把关:检验这些参数是否显著不为0 • 模型总体检验:
•检验1:𝐻:𝛽 =0 01
• 单个参数检验: •检验2:𝐻:𝛽 =0
00

模型回归系数评估(检验的角度)
መመ
• 问3:𝛽 与𝛽 服从什么分布?
01
• 答3:正态分布(前提是𝜀𝑖服从正态分布)
መመ
• 问4:给定检验2，3的原假设，𝛽 /𝑠 与𝛽 /𝑠 服从什么分布?
• 答4:自由度为𝑛 − 2的𝑡分布
0011

模型回归系数评估(检验的角度)
•检验1:𝐻:𝛽 =0 01
• H 成立时，模型退化为𝑦 = 𝛽 + 𝜖 0 𝑖0𝑖
• 此时其实一点信号都没有，如果我们当做有信号的话，就会有问题。

模型回归系数评估(检验的角度)
•
•
检验2:𝐻:𝛽 =0 00
检验3:𝐻:𝛽 =0 01
መመ 𝑖01
• 当𝜀 为正态分布时，𝛽 与𝛽 也是正态分布
෡
መ 𝛽 −𝛽
•𝛽∼𝑁𝛽,𝜎2 ⇔𝑗 𝑗∼𝑁(0,1)
𝜎𝑗
• 因为真实的𝜎 未知，只能用估计值𝑠 来代替 𝑗𝑗
෡
𝛽 −𝛽
• 检验3检验的是𝑥与𝑦之间的关系，检验2检验的是常数项是否有必要
𝑗𝑗𝑗
• 此时 𝑗 𝑗 ∼ 𝑡(𝑛 − 2)，如果这个值的绝对值特别大，则拒绝原假设 𝑠𝑗

多元线性回归 数据与模型

多元线性回归:数据
• 例:Country Kitchen公司零食年销售额
年销售额 广告投放额 促销投放额 竞争者年销售额 (百万 地区 (百万美元) (百万美元) (百万美元) 美元)
Selkirk 101.8 1.3 0.2 20.4
Susquehanna 44.4 0.7 0.2 30.5
Kittery 108.3 1.4 0.3 24.6
Acton 85.1 0.5 0.4 19.6
Finger Lakes 77.1 0.5 0.6 25.5
Berkshire 158.7 1.9 0.4 21.7
Central 180.4 1.2 1.0 6.8
Providence 64.2 0.4 0.4 12.6
Nashua 74.6 0.6 0.5 31.3
Dunster 143.4 1.3 0.6 18.6
Endicott 120.6 1.6 0.8 19.9
Five-Towns 69.7 1.0 0.3 25.6
Waldeboro 67.8 0.8 0.2 27.4
Jackson 106.7 0.6 0.5 24.3
Stowe 119.6 1.1 0.3 13.7
𝑌 𝑋 𝑋2 𝑋3 1

多元线性回归:数据
•模型:𝑦=𝛽+𝛽𝑥 +⋯+𝛽𝑥 +𝜖 𝑖=1,⋯,𝑛 𝑖01𝑖1 𝑝𝑖𝑝𝑖
• 𝜖𝑖 ∼𝑖.𝑖.𝑑. 𝑁(0, 𝜎2)
• 下图为𝑝 = 2的情况

多元线性回归:参数估计
•真实的模型是𝑦=𝛽+𝛽𝑥 +⋯+𝛽𝑥 +𝜖 𝑖01𝑖1 𝑝𝑖𝑝𝑖
መመመ •我们想找一条线去拟合:𝑦ො=𝛽+𝛽𝑥 +⋯+𝛽𝑥
መመ𝑛2 •最小二乘法:找𝛽 ,⋯,𝛽 ，使得σ (𝑦 −𝑦ො ) 最小
𝑖01𝑖1 𝑝𝑖𝑝
0𝑝 𝑖=1𝑖𝑖

多元线性回归:参数估计
CountryKitchen
= read.csv(“CountryKitchen.csv”,header=T,stringsAsFactors=F)
y = CountryKitchen
[,2]
x = as.matrix
(CountryKitchen[,-c(1,2)])
result2
= lm(y~x)
Sales =
CountryKitchen[,2]
Advertisment
= CountryKitchen[,3]
Promotion
= CountryKitchen[,4]
Competitor
= CountryKitchen[,5]
result3
= lm(Sales~Advertisment+Promotion+Competitor)
summary
(result3)

多元线性回归:参数估计
𝑌=𝛽+𝛽𝑥 +𝛽𝑥 +𝛽𝑥 +𝜖, 𝑖=1,⋯,𝑛 𝜖∼ 𝑁(0,𝜎2) 𝑖 0 1 1𝑖 2 2𝑖 3 3𝑖 𝑖 𝑖 𝑖.𝑖.𝑑.

多元线性回归:参数估计
判定系数𝑅2
𝑠 = 17.6

多元线性回归:参数估计
መ
𝛽 的估计值:𝛽 = 65.70
00
መ
𝛽 的估计值:𝛽 = 48.98
11
መ
𝛽 的估计值:𝛽 = 59.65
22
当所有的自变量为0时，公司年销售额期望的估 计值为65.70 百万美元
给定所有其他自变量保持不变时，每增加一个 单位 (1百万美元) 的广告投放额，公司期望增
加的年销售额估计值为48.98百万美元 给定所有其他自变量保持不变时，每增加一个
单位 (1百万美元) 的促销投放额，公司期望增 给定所有其他自变量保持不变时，竞争者年销
መ
𝛽 的估计值:𝛽 = −1.84
加的年销售额估计值为59.65百万美元 售额每增加一个单位 (1百万美元)，公司期望减
33
少的年销售额估计值为1.84百万美元

多元线性回归:参数估计
መ
𝛽 点估计
෢2 𝑠𝑖= 𝜎
假设检验
𝐻:𝛽=0 0𝑖
𝑖
መ
𝛽 的样本标准差
𝑖
𝑖
෡
𝛽𝑖−𝛽𝑖服从自由度为𝑛 − 𝑝 − 1的𝑡分布 𝑠𝑖

多元线性回归:模型的显著性检验
• 模型的显著性检验
•𝐻:𝛽=⋯=𝛽=0 v.s. 𝐻:至少有一个𝛽≠0
01𝑝𝑎𝑖
• 如果原假设成立，则表明模型中所有的p个自变量与应变量的取值均无关，也就是 说整个回归模型是不显著有效的
• p-值为 0.0001388 < 0.05，在0.05水平下拒绝 H0 模型参数评估(检验的角度) • 问:为什么下列检验不合在一起做? • 比如我检验了每一个𝛽 = 0，不就等于检验了𝛽 = ⋯ = 𝛽 = 0么? 𝑖1𝑝 多元线性回归 防忽悠进阶知识 吾日三省吾身 ◼𝑦=𝛽+𝛽𝑥+⋯+𝛽𝑥+𝜀 𝑖01𝑖1 𝑝𝑖𝑝𝑖 ◼ 模型评估问题(单一标准会产生的误区): ◼ 𝑅2是一个衡量模型拟合整体质量的完美度量吗? ◼ 模型(变量)选择问题: ◼ 如果遗漏了一些因素(变量𝑥)会怎么样? ◼ 我们多考虑一些因素(变量𝑥)会怎么样? ◼ 如果上面两个同时发生会怎么样? ◼ 模型诊断问题: ◼ 使用的模型是否与真实的情况吻合(适用性) 多元线性回归 关于𝑹𝟐进一步讨论 关于𝑹𝟐的进一步讨论 • 𝑅2是一个衡量模型拟合整体质量的完美度量吗? • Anscombe's Quartet 关于𝑹𝟐的进一步讨论 • 𝑅2是一个衡量模型拟合整体质量的完美度量吗? • 在数据中增加一个不相关的自变量: • 𝑥4:当年平均降雪量(英寸) 地区 年销售额 广告投放额 促销投放额 (百万美元) (百万美元) (百万美元) 竞争者年 年平均降雪量 销售额 (百 万美元) (英寸) Selkirk 101.8 1.3 0.2 20.4 24 Susquehanna 44.4 0.7 0.2 30.5 31 Kittery 108.3 1.4 0.3 24.6 31 Acton 85.1 0.5 0.4 19.6 36 Finger Lakes 77.1 0.5 0.6 25.5 18 Berkshire 158.7 1.9 0.4 21.7 42 Central 180.4 1.2 1 6.8 50 Providence 64.2 0.4 0.4 12.6 49 Nashua 74.6 0.6 0.5 31.3 60 Dunster 143.4 1.3 0.6 18.6 62 Endicott 120.6 1.6 0.8 19.9 42 Five-Towns 69.7 1 0.3 25.6 58 Waldeboro 67.8 0.8 0.2 27.4 55 Jackson 106.7 0.6 0.5 24.3 79 Stowe 119.6 1.1 0.3 13.7 88 关于𝑹𝟐的进一步讨论 Snow = c(24,31,31,36,18,42,50,49,60,62,42,58,55,79,88) result4 = lm(Sales~Advertisment+Promotion+Competitor+Snow) summary(result4) 关于𝑹𝟐的进一步讨论 •𝑌=𝛽+𝛽𝑥 +𝛽𝑥 +𝛽𝑥 +𝜀 𝑖 0 11𝑖 22𝑖 33𝑖 𝑖 •𝑌=𝛽+𝛽𝑥 +𝛽𝑥 +𝛽𝑥 +𝛽𝑥 +𝜀 𝑖 0 11𝑖 22𝑖 33𝑖 44𝑖 𝑖 • 𝑅2从0.833增加到0.856 • 基于𝑅2，第二个模型似乎比第一个模型更好，但是从直观上来看𝑅2与年销售额无关，并不应该将其 放入模型中 • 一个事实:𝑅2总是随着新的自变量不断进入模型而增加，哪怕该自变量与应变量几乎没有关系，𝑅2 也会小幅增加 • 然而在建模的过程中，不应该只为了使得𝑅2在表象上增加而将与应变量无关的自变量引入模型 关于𝑹𝟐的进一步讨论 ⚫𝑅2的一个问题: ⚫𝑅2总是随着自变量个数𝑝 的增加而增加 ⚫ 在一个模型中引入多个与应变量无关的自变量时也可以在表面上提 高𝑅2 ⚫ 统计学家们提出了不少方法来克服判定系数的这一缺陷，其中 一种是使用修正判定系数来衡量模型的整体拟合质量 𝑅2 =1−(1−𝑅2)( 𝑛−1 ) 𝑎 𝑛−𝑝−1 ⚫ 相比于判定系数本身，修正判定系数还考虑了模型中使用的自 变量个数𝑝与样本容量𝑛 关于𝑹𝟐的进一步讨论 •𝑌=𝛽 +𝛽𝑥 +𝛽𝑥 +𝛽𝑥 +𝜀 𝑖 0 11𝑖 22𝑖 33𝑖 𝑖 •𝑌=𝛽 +𝛽𝑥 +𝛽𝑥 +𝛽𝑥 +𝛽𝑥 +𝜀 𝑖 0 11𝑖 22𝑖 33𝑖 44𝑖 𝑖 • 𝑅2的增加量比𝑅2的增加量小，这是因为前者考虑了模型中自变量x的个数 𝑎 • 然而在这个例子中，即便是基于𝑅2进行判断，仍会将无应变量几乎无关的自变量引 𝑎 入模型，修正判定系数只是一种修正的方法，未必在任何情况下都有效 • 为了衡量模型中每个自变量是否确实与应变量有关，一种更正规和有效的方法是对 回归系数的估计进行显著性检验 关于𝑹𝟐的进一步讨论 • 加入无关的自变量𝑥4:年平均降雪量 • 在0.05的水平下拒绝 𝐻 :𝛽 =𝛽 =𝛽 =𝛽 =0 01234 • 在0.05的水平下不能拒绝 𝐻 : 𝛽 = 0 04 • 显著性检验表明年平均降雪量这一自变量对应变量取值没有统计意义下的显著影响， 因此应该从模型中除去 多元线性回归 多重共线性与变量选择 多重共线性 ◼ 例:学校表现与起薪 学生编号 起薪 ($) 大学GPA GMAT 1 100,000 3.9 640 2 100,000 3.9 644 3 77,500 3.1 557 4 77,500 3.2 550 5 75,000 3.0 547 6 87,500 3.5 589 7 77,500 3.0 533 8 87,500 3.5 600 9 77,500 3.2 630 10 90,000 3.6 633 11 95,000 3.7 642 12 65,000 2.9 522 13 72,500 3.4 628 14 82,500 3.4 583 15 100,000 4.0 650 16 97,500 3.8 641 17 68,500 2.8 530 18 85,000 3.5 596 19 88,500 3.6 605 20 100,000 3.9 656 21 77,500 3.1 574 22 92,500 3.7 636 23 92,500 3.7 635 24 97,500 4.0 654 25 95,000 3.8 633 Y:起薪 x1 : 大学GPA x2 : GMAT 多重共线性 (, ,,,,,,, ,,,,,, salary=c 100000 95000 , , 65000, 100000, 72500, 77500, 82500, 77500 100000 75000 97500 87500 68500 77500 85000 87500 88500 77500 100000 90000 77500 92500, 92500, 97500, 95000 GPA=c ( , , , , 3.9 , , , , , 3.9 , , , , , 3.1, 3.2 3.0 3.5 3.0 3.5 3.2, 3.6 3.7 2.9 3.4 3.4, 4.0 3.8 2.8 3.5 3.6, 3.9 3.1 3.7 3.7 4.0, 3.8 GMAT=c ( 640 , 644, 557 550 589 , , , , 533 , , , , ( 600 ,6 ,5 ,6 ,6 30 , , , , 633 642 522 628 83 650 641 530 596 05 656 574 636 635 54 633 result5 = lm salary~GPA+GMAT summary result5 ) ,, , , , ) , ,547, , , , ) ) () 多重共线性 • 判定系数与修正判定系数表示模型拟合度不错 • GMAT成绩对应的系数估计值与预期不一致，而且对应的𝑝 − 𝑣𝑎𝑙𝑢𝑒显示 其不显著。 • 一定是哪里出了问题! 多重共线性 • 看一下三个变量的样本相关系数 起薪 大学GPA GMAT 起薪 1 大学GPA 0.955642 1 GMAT 0.847114 0.915028 1 ⚫ 两个自变量有很高的相关系数，表明两者之间按有明显的线性关系 ⚫ 这样的两个自变量同时出现在一个回归模型中，它们的作用可以彼此替代 ⚫ 这是典型的多重共线性 ⚫ 多重共线性问题常常可以通过从回归模型中删除一个或多个自变量得以缓解 ⚫ 一般会删除与因变量相关性最弱的自变量 多重共线性 • 从模型中删除GMAT 变量选择 ⚫ 现在已可以看出在回归模型中引入过多自变量可能会引起严重的问题， 因此在所有可能的自变量中，应只选择一些与应变量最“有关”的进行 建模，这就是“变量选择”问题 ⚫ 变量选择可以缓解多重共线性 变量选择 例:Cravens数据 ⚫ Cravens数据是一家公司的数据，这家公司在一些销售区域销售产品，并且为每 一个销售区域指定了独家经销商，回归分析的目的是为了确定:预测(自)变量 的变化能否解释每一个区域的销售情况;由25个销售区域组成了一个随机样本， 得到数据: 变量选择 ⚫初步分析:变量之间的样本相关系数 Sales Time Poten AdvExp Share Change Accounts Work Rating Sales 1 Time Poten 0.622881 1 0.597812 0.453978 1 AdvExp 0.596177 0.24917 0.1741 1 Share 0.483511 0.106105 − 0.21067 0.264458 1 Change 0.489181 0.251469 0.268287 0.376516 0.08547 1 Accounts Work 0.753986 0.757789 0.478643 0.200037 0.403011 0.327415 1 − 0.11722 − 0.17941 − 0.25884 − 0.27223 0.349345 − 0.28766 − 0.19885 1 Rating 0.401878 0.101172 0.358703 0.411462 − 0.02356 0.549333 0.228613 − 0.27691 1 ◼哪些自变量与应变量有较强的样本相关性?—— 应该选择哪些自变量来 预测 ◼ 哪些自变量之间有较强的样本相关性?——有没有潜在的多重共线性 ◼ 变量选择:一种选择自变量的自动算法 变量选择 ◼ 可以用p-值作为选择变量的准则——剔除大于0.05的自变量 这些自变量 似乎不应留 在模型中? 至此根据p-值已不能再剔除自变量，从8个自变量中选择了4个，剩 余的4个自变量在0.05水平下均显著 变量选择 ◼ 变量选择可以缓解多重共线性 ◼ 例:Cravens数据 自变量Time与Accounts有较高的样本相关系数(0.76)，有存在多重共线性的危险 两个自变量都与应变量 有较高的样本相关系数， 但同时出现回归模型中 时p-值都大于0.05，这 是多重共线性的特征之 一 变量选择有效地减轻多 重共线性带来的问题 变量选择 ⚫ 一般地，上述的变量选择过程可以总结如下: ⚫ 假设总共有 p 个可能的自变量，用所有的 p 个自变量拟合一个多元回归模型 ; ⚫ 如果其中有一些自变量系数估计值的p-值大于0.05，则找出其中p-值最大的自 变量剔除出模型，用余下的 p − 1个自变量拟合多元回归模型; ⚫ 重复步骤 2 直至基于p-值没有自变量可以从回归模型中可以被剔除 ⚫ 上述的变量选择过程被称为向后消元法 ⚫ 还有另外两种常用的变量选择方法:向前选择法与逐步回归法 ⚫ 常用的统计软件中都有这三种变量选择方法 变量选择(R中实现) Cravens = read.csv("Cravens.csv",header=T,stringsAsFactors=F) Sales Time Poten AdvExp = Cravens[,4] Share = Cravens[,5] Change = Cravens[,6] Accounts= Cravens[,7] Work = Cravens[,8] Rating = Cravens[,9] lm7 =lm(Sales~Time+Poten+AdvExp+Share+Change+Accounts+Work+Rating) step(lm7,direction="back") 注意:R中选择的标准为AIC = Cravens[,1] = Cravens[,2] = Cravens[,3] 变量选择(R中实现) 变量选择 ⚫ 如果分析的时候遗漏重要变量，会对结果产生严重的误判 变量选择 ⚫ 如果分析的时候遗漏重要变量，会对结果产生严重的误判 变量选择 ⚫ 如果分析的时候遗漏重要变量，会对结果产生严重的误判 多元线性回归 如何使用 使用回归模型:描述 ⚫ 回归模型的主要功能有两项:描述与预测 ⚫ 关于描述:我们基于回归系数的最小二乘估计值，给出自变量对于应变 量取值影响的方式与影响的大小 (估计值的解释) ⚫关于预测:现有一组新的自变量值 (𝑥∗,⋯,𝑥∗)，但其对应的应变量取值 𝑦∗为多少? 1𝑝 使用回归模型:预测 ⚫ 关于预测:可以基于拟合的回归方程对一个新的、应变量取值为观测的 个体进行预测 ⚫ 假设回归的模型为 ⚫ 拟合的回归方程为 ⚫ 假设此个体仍满足多元线性回归模型 𝑦∗=𝛽 +𝛽𝑥∗+⋯+𝛽𝑥∗+𝜀∗ ⚫ 对于因变量𝑦∗取值的点估计为: 𝑦=𝛽+𝛽𝑥+⋯+𝛽𝑥+𝜀 𝑖01𝑖1 𝑝𝑖𝑝𝑖 መመመ 𝑦ො = 𝛽 + 𝛽 𝑥 + ⋯ + 𝛽 𝑥 𝑖01𝑖1 𝑝𝑖𝑝 1𝑝 ⚫现有一组新的自变量值 (𝑥∗,⋯,𝑥∗)，但其对应的应变量取值 𝑦∗未知 011𝑝𝑝 ∗መመ∗መ∗ 𝑦ො = 𝛽 + 𝛽 𝑥 + ⋯ + 𝛽 𝑥 011𝑝𝑝 使用回归模型:预测 ⚫Country Kitchen公司零食年销售额 𝑦ො𝐶𝑜𝑙𝑜𝑚𝑏𝑢𝑠 =65.7+49.0×1.2+59.7×0.5−1.8×25=109.35 𝑦ො𝐵𝑙𝑜𝑜min𝑔𝑡𝑜𝑛 =65.7+49.0×1+59.7×0.6−1.8×10=132.52 𝑦ො𝐵𝑢𝑓𝑓𝑎𝑙𝑜 =65.7+49.0×1.4+59.7×0.7−1.8×8=161.69 地区 广告投放额 促销投放额 (百万美元) (百万美元) 竞争者年销售 年销售额 额 (百万美元) (百万美元) Columbus 1.2 0.5 25 109.35 Bloomington 1 0.6 10 132.52 Buffalo 1.4 0.7 8 161.69 多元线性回归 模型诊断 回归诊断 ⚫ 自圆其说很重要! ⚫ 以下这些假定在真实的数据中满足么? ⚫ 线性回归模型 ⚫𝑦=𝛽+𝛽𝑥+⋯+𝛽𝑥+𝜀 𝑖01𝑖1 𝑝𝑖𝑝𝑖 ⚫ 上述多元线性回归模型的假设包括 ⚫ 每个自变量对于应变量取值的影响均服从线性关系 ⚫ 随机误差与自变量独立 ⚫ 随机误差相互独立(独立性) ⚫ 随机误差的分布相同，均值为0，并且有相同的标准差 (同方差性) ⚫ 随机误差服从正态分布 回归诊断 ⚫ 在把拟合好的回归模型进行应用之前，还需要对拟合的回归模型进行诊 断，以判断那些对于模型所施加的假设是否与拟合结果一致 ⚫ 假设在进行了显著性检验与变量选择之后，在回归模型中最后保留了𝑝 个自变量来解释应变量，记估计的回归方程为 መመመ 𝑦ො = 𝛽 + 𝛽 𝑥 + ⋯ + 𝛽 𝑥 𝑖01𝑖1 𝑝𝑖𝑝 ⚫ 假定真实的模型为 𝑦=𝛽+𝛽𝑥+⋯+𝛽𝑥+𝜀 ⚫ 残 差 : 𝑒 = 𝑦 − 𝑦ො 𝑖𝑖𝑖 ⚫当回归成立时:𝑒 ,⋯,𝑒 是随机误差𝜀 ,⋯,𝜀 的估计 1𝑛 1𝑛 𝑖01𝑖1 𝑝𝑖𝑝𝑖 检测残差与自变量独立性 例:Country Kitchen公司零食年销售额 30 20 10 30 20 10 广告投放额 (百万美元) 残差图 竞争者年销售额 (百万美元) 残差图 0 -100 10 20 30 40 -20 -30 -40 竞争者年销售额 (百万美元) 0 -100 0.5 1 1.5 2 -20 -30 -40 30 20 10 0 广告投放额 (百万美元) 促销投放额 (百万美元) 残差图 0 0.2 0.4 0.6 0.8 1 1.2 -10 -20 -30 -40 促销投放额 (百万美元) 残差 残差 残差 检测异方差性 ⚫ 检测异方差性:画出残差(纵轴)对单个自变量或拟合应变量(横轴)的散 点图 (残差图)，从图中观察残差对于单个自变量或拟合应变量的取值 是否呈现出(函数形式的)趋势 56 4 3 2 1 5 4 3 2 1 00 -10 0.5 1 1.5 2-10 0.5 1 1.5 2 -2 -3 -4 -2 -3 -4 -5 拟合的应变量 无异方差性的证据 拟合的应变量 有异方差性的证据 残差 残差 检测残差之间独立性(最难检验) ◼ 自相关性是一种常见的违反独立性的现象，当样本观测是以某种自然序 列(如时间)被记录是，可能出现自相关性。 ◼ 将残差按自然序列作图并检测趋势 无自相关性证据 有自相关性证据 检测正态性 ⚫ 随机误差的正态性假设可以用残差的 ◼ 如果直方图基本上呈现出对称与倒钟型的形态，则表明没有明显证据违 背正态性假设 ◼ 然而，当样本容量𝑛很小时直方图作用较小 直方图来进行检测 检测正态性 ⚫例:Country Kitchen公司零食年销售额 Q-Q Normal Plot of Residuals 1.50 0.50 -2.00 -1.00 0.00 1.00 2.00 -0.50 -1.50 -2.50 ◼ 图中的点越靠近直线，则正态 性假设越充分 ◼ 根据由例中数据得到的Q-Q图， 没有明显的证据显示随机误差 不服从正态分布 标准化残差 多元线性回归 回归诊断案例 Courtesy of Prof. Xuming He 回归诊断案例 • Time Series Regression (Hotel Occupancy Case) • In a recent legal case, a Chicago downtown hotel claimed that it had suffered a loss of business due to what was considered an illegal action by a group of hotels that decided to leave the plaintiff out of a hotel directory. • In order to estimate the loss business, the hotel had to predict what its level of business (in terms of occupancy rate) would have been in the absence of the alleged illegal action. • In order to do this, experts testifying on behalf of the hotel use data collected before the period in question and fit a relationship between the hotels occupancy rate and overall occupancy rate in the city of Chicago. This relationship would then be used to predict occupancy rate during the period in question. Courtesy of Prof. Xuming He 回归诊断案例 回归诊断案例 回归诊断案例 回归诊断案例 回归诊断案例 回归诊断案例 回归诊断案例 回归诊断案例 回归诊断案例 含有定性自变量的回归模型 线性回归模型 •模型:𝑦=𝛽+𝛽𝑥 +⋯+𝛽𝑥 +𝜀 𝑖01𝑖1 𝑝𝑖𝑝𝑖 • 𝑖 = 1, ⋯ , 𝑛 • 所有的自变量与应变量都是连续型变量 • 如果自变量与应变量其中一个，或者全部，都是离散型的怎么办? • 比如: • 如果 𝑥 是0-1型的，或是多类/定序的(好、中、差) • 如果 𝑦 是0-1型的，或是多类/定序的(好、中、差) • 如果 𝑦 是一个计数类型的 𝑿为定性变量， 𝒀为连续变量 自变量为定性变量 ◼ 在许多情形下，自变量为定性变量 (名义型或有序型) ◼ 例 :约翰逊过滤水股份公司 约翰逊过滤水股份有限公司对于遍布南弗罗里达州的水过滤系统进行维修;为 了估计服务时间和服务成本，公司的管理人员希望对顾客的每一次维修请求预测必 要的维修时间，管理人员认为维修时间依赖于两个因素:从最近一次维修服务至今 水过滤系统已经使用的月数和需要维修的故障类型 (机械的或电子的) 维修服务请求 维修时间 上次维修至今的月数 维修的故障类型 (小时) 1 2.9 2 电子 2 3.0 6 机械 3 4.8 8 电子 4 1.8 3 机械 5 2.9 2 电子 6 4.9 7 电子 7 4.2 9 机械 8 4.8 8 机械 9 4.4 4 电子 10 4.5 6 电子 虚拟变量法 ◼ 自变量“维修的故障类型”是一个名义型变量 ◼ 如何将其放入回归模型中? ◼ 引入一个“虚拟变量”来 量化该自变量 维修服务请求 维修时间 上次维修至今的月 维修的故 (小时) 数 障类型 故障类型 1 2.9 2 电子 1 2 3.0 6 机械 0 3 4.8 8 电子 1 4 1.8 3 机械 0 5 2.9 2 电子 1 6 4.9 7 电子 1 7 4.2 9 机械 0 8 4.8 8 机械 0 9 4.4 4 电子 1 10 4.5 6 电子 1 虚拟变量法 hours = c(2.9, 3.0, 4.8, 1.8, 2.9, 4.9, 4.2, 4.8, 4.4, 4.5) gap = c(2, 6, 8, 3, 2, 7, 9, 8, 4, 6) error_type = c(1, 0, 1, 0, 1, 1, 0, 0, 1, 1) lm1 = lm(hours ~ 1+gap+error_type) summary(lm1) 如何解读? 虚拟变量法 𝑦ො = 0.93 + 0.39𝑔𝑎𝑝 + 1.26 × 𝑒𝑟𝑟_𝑡𝑦𝑝𝑒 ◼ 机械故障:𝑦ො = 0.93 + 0.39𝑔𝑎𝑝 ◼ 电子故障:𝑦ො = 0.93 + 0.39𝑔𝑎𝑝 + 1.26 ◼ 给定其他自变量不变的情况下，电子故障的维修时间要比机械 故障的维修时间平均要多1.26个小时。 𝑿为多类型定性变量， 𝒀为连续变量 金融投资收益率的分析 金融投资收益率的分析 ◼ 金融投资有很多种，不同的投资对应着不同的收益率(期望)，当然， 也对应着不同的风险(方差)。 ◼ 我们暂时只考虑收益率。 ◼ 有些投资设置了一定的门槛，只有当资金超过一定的量才允许参加投 资。 ◼ 这些投资是不同类型的银行运营的，如国有控股银行，股份制商业银 行，城市商业银行，其他类型银行。 ◼ 我们感兴趣的是，不同投资门槛，不同类型银行所运营的投资在收益 率上有什么不同。 金融投资收益率的分析 ◼ 部分数据 金融投资收益率的分析 ◼ 我们对数据重新编码: ◼ 投资门槛(以低门槛为基准) ◼𝑋 =0 (低门槛) 1 ◼𝑋 =1 (高门槛) 1 ◼ 银行类型(以国有控股银行为基准) ◼ 𝑋2 = 0 (非股份制商业银行) ◼ 𝑋2 = 1 (股份制商业银行) ◼ 𝑋3 = 0 (非城市商业银行) ◼ 𝑋3 = 1 (城市商业银行) ◼ 𝑋4 = 0 (非其他类型银行) ◼ 𝑋4 = 1 (其他类型银行) 金融投资收益率的分析 ◼ 编码后部分数据 金融投资收益率的分析 ◼模型:𝑦 =𝛽 +𝛽𝑥 +𝛽𝑥 +𝛽𝑥 +𝛽𝑥 +𝜀 𝑖011223344𝑖 ◼ 如何建立模型? file_name = "net_value.csv" data1 = read.csv(file = file_name) ret = data1[,3] qidian = data1[,4]; joint_equity_bank = data1[,5]; city_commerical_bank = data1[,6]; other_bank = data1[,7] lm1 = lm(ret ~ 1+qidian+joint_equity_bank+city_commerical_bank+other_bank) summary(lm1) 金融投资收益率的分析 ◼模型:𝑦 =𝛽 +𝛽𝑥 +𝛽𝑥 +𝛽𝑥 +𝛽𝑥 +𝜀 𝑖011223344𝑖 ◼ 如何解读模型分析结果? 金融投资收益率的分析 ◼ 投资门槛(以低门槛为基准) ◼𝑋 =0 (低门槛) ◼𝑋 =1 (高门槛) 1 ◼ 银行类型(以国有控股银行为基准) ◼ 𝑋2 = 0 (非股份制商业银行) ◼ 𝑋2 = 1 (股份制商业银行) ◼ 𝑋3 = 0 (非城市商业银行) ◼ 𝑋3 = 1 (城市商业银行) ◼ 𝑋4 = 0 (非其他类型银行) ◼ 𝑋4 = 1 (其他类型银行) 𝑦=𝛽+𝛽𝑥+𝛽𝑥+𝛽𝑥+𝛽𝑥+𝜀 𝑖011223344𝑖 1 低门槛 高门槛 国有控股银行 𝑦=𝛽+𝜀 𝑖0𝑖 𝑦=𝛽+𝛽+𝜀 𝑖01𝑖 股份制银行 𝑦=𝛽+𝛽+𝜀 𝑖02𝑖 𝑦=𝛽+𝛽+𝛽+𝜀 𝑖012𝑖 城市商业银行 𝑦=𝛽+𝛽+𝜀 𝑖03𝑖 𝑦=𝛽+𝛽+𝛽+𝜀 𝑖013𝑖 其他类型银行 𝑦=𝛽+𝛽+𝜀 𝑖04𝑖 𝑦=𝛽+𝛽+𝛽+𝜀 𝑖014𝑖 金融投资收益率的分析 ◼ 我们对数据重新编码: ◼ 投资门槛(以低门槛为基准) ◼𝑋 =0 (低门槛) 1 ◼𝑋 =1 (高门槛) 1 ◼ 银行类型(以国有控股银行为基准) ◼ 𝑋2 = 0 (非股份制商业银行) ◼ 𝑋2 = 1 (股份制商业银行) ◼ 𝑋3 = 0 (非城市商业银行) ◼ 𝑋3 = 1 (城市商业银行) ◼ 𝑋4 = 0 (非其他类型银行) ◼ 𝑋4 = 1 (其他类型银行) ◼ 𝑋5 = 0 (非国有控股银行) ◼ 𝑋5 = 1 (国有控股银行) 如果国有控股银行也作为 一个变量呢? 金融投资收益率的分析 ◼𝑦=𝛽+𝛽𝑥+𝛽𝑥+𝛽𝑥+𝛽𝑥+𝛽𝑥+𝜀 𝑖01122334455𝑖 ◼ 且𝑥2 + 𝑥3 + 𝑥4 + 𝑥5 = 1(𝑥2, 𝑥3, 𝑥4, 𝑥5同时有且只有一个为1，其他的都为0) ◼𝑦=𝛽+𝛽𝑥+𝛽𝑥+𝛽𝑥+𝛽𝑥+𝛽1−𝑥−𝑥−𝑥 +𝜀 𝑖0112233445 234𝑖 ◼𝑦=𝛽+𝛽+𝛽𝑥+𝛽−𝛽 𝑥+𝛽−𝛽 𝑥+(𝛽−𝛽)𝑥+𝜀 𝑖0511252353454𝑖 讨论:金融投资收益率的分析 𝑦=𝛽+𝛽𝑥+𝛽𝑥+𝜀 𝑖01122𝑖 ◼ 银行类型 ◼ 𝑥2 = 0 (国有控股银行) ◼ 𝑥2 = 1 (股份制商业银行) ◼ 𝑥2 = 2 (城市商业银行) ◼ 𝑥2 = 3 (其他类型银行) ◼ 银行类型 ◼国有控股银行:𝑦 =𝛽 +𝛽𝑥 +𝜀 𝑖011𝑖 ◼股份制商业银行:𝑦 =𝛽 +𝛽𝑥 +𝛽 +𝜀 𝑖0112𝑖 ◼城市商业银行:𝑦 =𝛽 +𝛽𝑥 +2𝛽 +𝜀 𝑖0112𝑖 ◼其他类型银行:𝑦 =𝛽 +𝛽𝑥 +3𝛽 +𝜀 𝑖0112𝑖 𝑿为定序变量， 𝒀为连续变量 北京市某年房价与其位置关系 ◼ 部分数据 位置 价格 四至五环 6200 四至五环 6200 四至五环 6300 四至五环 6480 四至五环 7200 四至五环 7800 四至五环 7800 四至五环 6100 四至五环 3800 四至五环 5880 四至五环 7800 四至五环 4756 四至五环 5800 四至五环 6200 北京市某年房价与其位置关系 ◼ 我们对数据重新编码: ◼𝑌 价格 ◼ 位置 ◼𝑋 是否在二环以内?(是则为1，否则为0) 1 ◼𝑋2 是否在二至三环?(是则为1，否则为0) ◼𝑋3 是否在三至四环?(是则为1，否则为0) ◼𝑋4 是否在四至五环?(是则为1，否则为0) ◼𝑋5 是否在五环以外?(是则为1，否则为0) 北京市某年房价与其位置关系 ring X1 X2 X3 X4 X5 price 四至五环 0 0 0 1 0 6200 四至五环 0 0 0 1 0 6200 四至五环 0 0 0 1 0 6300 四至五环 0 0 0 1 0 6480 四至五环 0 0 0 1 0 7200 四至五环 0 0 0 1 0 7800 四至五环 0 0 0 1 0 7800 四至五环 0 0 0 1 0 6100 四至五环 0 0 0 1 0 3800 四至五环 0 0 0 1 0 5880 四至五环 0 0 0 1 0 7800 四至五环 0 0 0 1 0 4756 四至五环 0 0 0 1 0 5800 四至五环 0 0 0 1 0 6200 北京市某年房价与其位置关系 ◼ 我们对数据重新编码: ◼ 𝑌 价格 ◼ 位置 ◼𝑋 是否在二环以内?(是则为1，否则为0) 1 ◼𝑋2 是否在二至三环?(是则为1，否则为0) ◼𝑋3 是否在三至四环?(是则为1，否则为0) ◼𝑋4 是否在四至五环?(是则为1，否则为0) ◼𝑋5 是否在五环以外?(是则为1，否则为0) ◼ 如何建模? ◼𝑌=𝛽 +𝛽𝑋 +𝛽𝑋 +𝛽𝑋 +𝛽𝑋 +𝜀 011223344 ◼𝑌=𝛽 +𝛽𝑋 +𝛽𝑋 +𝛽𝑋 +𝛽𝑋 +𝛽𝑋 +𝜀 01122334455 北京市某年房价与其位置关系 file_name = "BJ_REAL_1.csv" data1 = read.csv(file = file_name, stringsAsFactors = F) ring = data1[,1] X1 = data1[,2] X2 = data1[,3] X3 = data1[,4] X4 = data1[,5] X5 = data1[,6] price = data1[,7] lm1 = lm(price ~ 1+X1+X2+X3+X4) summary(lm1) lm2 = lm(price ~ 1+X1+X2+X3+X4+X5) summary(lm2) 北京市某年房价与其位置关系 北京市某年房价与其位置关系 𝑌=𝛽 +𝛽𝑋 +𝛽𝑋 +𝛽𝑋 +𝛽𝑋 +𝜀 011223344 ◼ 位置 ◼𝑋 是否在二环以内?(是则为1，否则为0) 1 ◼𝑋2 是否在二至三环?(是则为1，否则为0) ◼𝑋3 是否在三至四环?(是则为1，否则为0) ◼𝑋4 是否在四至五环?(是则为1，否则为0) ◼𝑋5 是否在五环以外?(是则为1，否则为0) ◼ 如果房子在 ◼二环以内:𝑌=𝛽 +𝛽 +𝜀 ◼二至三环:𝑌=𝛽 +𝛽 +𝜀 02 ◼三至四环:𝑌=𝛽 +𝛽 +𝜀 03 ◼四至五环:𝑌=𝛽 +𝛽 +𝜀 04 ◼五环以外:𝑌=𝛽 +𝜀 0 01 北京市某年房价与其位置关系 ◼ 如何解释? ◼ 𝑌 = 5074+5065𝑋 +4897𝑋 +3630𝑋 +1549𝑋 +𝜀 ◼ 五环之外均价5074元/平米 ◼ 四至五环比五环之外高1549元/平米 ◼ 三至四环比五环之外高3630元/平米 ◼ 二至三环比五环之外高4897元/平米 ◼ 二环以内比五环之外高5065元/平米 1234 北京市某年房价与其位置关系 ◼ 如何解释? ◼ 𝑌 = 5074+5065𝑋 +4897𝑋 +3630𝑋 +1549𝑋 +𝜀 1234 ◼ 五环之外均价5074元/平米 ◼ 四至五环比五环之外高1549元/平米 ◼ 三至四环比五环之外高3630元/平米 ◼ 二至三环比五环之外高4897元/平米 ◼ 二环以内比五环之外高5065元/平米 ◼ 改进: ◼ 为什么不考虑四至五环比五环之外高多少? ◼ 为什么不考虑三至四环比四至五环高多少? ◼ 为什么不考虑二至三环比三至四环高多少? ◼ 为什么不考虑二环以内比二至三环高多少? 北京市某年房价与其位置关系 ◼ 我们对数据重新编码: ◼𝑌 价格 ◼ 位置 ◼𝑋 是否在二环以内?(是则为1，否则为0) 1 ◼𝑋2 是否在三环以内?(是则为1，否则为0) ◼𝑋3 是否在四环以内?(是则为1，否则为0) ◼𝑋4 是否在五环以内?(是则为1，否则为0) ◼𝑋5 是否在五环以外?(是则为1，否则为0) 北京市某年房价与其位置关系 ring 二环以内 三环以内 四环以内 五环以内 五环之外 price 四至五环 0 0 0 1 0 6200 四至五环 0 0 0 1 0 6200 四至五环 0 0 0 1 0 6300 四至五环 0 0 0 1 0 6480 四至五环 0 0 0 1 0 7200 四至五环 0 0 0 1 0 8000 四至五环 0 0 0 1 0 13000 四至五环 0 0 0 1 0 8500 三至四环 0 0 1 1 0 7300 三至四环 0 0 1 1 0 18000 三至四环 0 0 1 1 0 10480 三至四环 0 0 1 1 0 8300 二至三环 0 1 1 1 0 10000 二至三环 0 1 1 1 0 12000 北京市某年房价与其位置关系 ◼ 我们对数据重新编码: ◼𝑌 价格 ◼ 位置 ◼𝑋 是否在二环以内?(是则为1，否则为0) 1 ◼𝑋2 是否在三环以内?(是则为1，否则为0) ◼𝑋3 是否在四环以内?(是则为1，否则为0) ◼𝑋4 是否在五环以内?(是则为1，否则为0) ◼𝑋5 是否在五环以外?(是则为1，否则为0) ◼ 如何建模? ◼𝑌=𝛽 +𝛽𝑋 +𝛽𝑋 +𝛽𝑋 +𝛽𝑋 +𝜀 011223344 ◼𝑌=𝛽 +𝛽𝑋 +𝛽𝑋 +𝛽𝑋 +𝛽𝑋 +𝛽𝑋 +𝜀 01122334455 北京市某年房价与其位置关系 北京市某年房价与其位置关系 ◼ 位置 ◼𝑋 是否在二环以内?(是则为1，否则为0) 1 ◼𝑋2 是否在三环以内?(是则为1，否则为0) ◼𝑋3 是否在四环以内?(是则为1，否则为0) ◼𝑋4 是否在五环以内?(是则为1，否则为0) ◼𝑋5 是否在五环以外?(是则为1，否则为0) ◼ 如果房子在 ◼ 二环以内:𝑌=𝛽 +𝛽 +𝛽 +𝛽 +𝛽 +𝜀 01234 ◼ 二至三环:𝑌=𝛽 +𝛽 +𝛽 +𝛽 +𝜀 0234 ◼三至四环:𝑌=𝛽 +𝛽 +𝛽 +𝜀 034 ◼四至五环:𝑌=𝛽 +𝛽 +𝜀 04 ◼五环以外:𝑌=𝛽 +𝜀 0 北京市某年房价与其位置关系 ◼ 位置 ◼𝑋 是否在二环以内?(是则为1，否则为0) 1 ◼𝑋2 是否在三环以内?(是则为1，否则为0) ◼𝑋3 是否在四环以内?(是则为1，否则为0) ◼𝑋4 是否在五环以内?(是则为1，否则为0) ◼𝑋5 是否在五环以外?(是则为1，否则为0) ◼ 如果房子在 ◼ 二环以内:𝑌=𝛽 +𝛽 +𝛽 +𝛽 +𝛽 +𝜀 01234 ◼二至三环:𝑌=𝛽 + 𝛽 +𝛽 +𝛽 +𝜀 0 234 ◼三至四环:𝑌=𝛽 + 𝛽 +𝛽 +𝜀 034 ◼四至五环:𝑌=𝛽 + 𝛽 +𝜀 04 ◼ 五环以外:𝑌 = 𝛽 + 𝜀 0 北京市某年房价与其位置关系 ◼ 如何解释? ◼ 𝑌 = 5074+168𝑋 +1267𝑋 +2081𝑋 +1549𝑋 +𝜀 ◼ 五环之外均价5074元/平米 ◼ 四至五环比五环之外高1549元/平米 ◼ 三至四环比四至五环高2081元/平米 ◼ 二至三环比三至四环高1267元/平米 ◼ 二环之内比二至三环高168元/平米 1234 北京市某年房价与其位置关系 ◼ 对比一下(系数): ◼ 五环之外均价5074元/平米 ◼ 四至五环比五环之外高1549元/平米 ◼ 三至四环比五环之外高3630元/平米 ◼ 二至三环比五环之外高4897元/平米 ◼ 三至四环比五环之外高5065元/平米 ◼ 五环之外均价5074元/平米 ◼ 四至五环比五环之外高1549元/平米 ◼ 三至四环比四至五环高2081元/平米 ◼ 二至三环比三至四环高1267元/平米 ◼ 二环之内比二至三环高168元/平米 北京市某年房价与其位置关系 ◼ 对比一下(显著性)，说明什么问题? 内环与外环比 都与五环外比 北京市某年房价与其位置关系 ◼ 看看地图 补充 交叉效应 ◼ 北京的空气质量，有人说是因为供暖，有人说是因为汽车排气。 ◼ 如果两件事情凑在一起怎么办? ◼ 定义: ◼ 𝑥1:供暖与否(定性变量) ◼ 𝑥2:出行车量(连续变量) ◼𝑌=𝛽+𝛽𝑥+𝛽𝑥+𝛽𝑥𝑥+⋯+ε 0 11 22 1212 对数线性回归模型 •𝑦=𝛽+𝛽𝑥+𝜀 𝑖01𝑖𝑖 • 𝑑𝑦 = 𝛽 𝑑𝑥 𝑖1𝑖 •𝑦=𝛽+𝛽log𝑥 +𝜀 𝑖01𝑖𝑖 • 𝑑𝑦 = 𝛽 𝑑𝑥𝑖 𝑖 1𝑥𝑖 • log(𝑦 ) = 𝛽 + 𝛽 𝑥 + 𝜀 𝑑𝑦 𝑖 0 1𝑖 𝑖 • 𝑖=𝛽𝑑𝑥 𝑦𝑖 1 𝑖 •log(𝑦)=𝛽 +𝛽log(𝑥)+𝜀 𝑖01𝑖𝑖 • 𝑑𝑦𝑖 = 𝛽 𝑑𝑥𝑖 𝑦𝑖 1 𝑥𝑖 回归分析小结 回归分析的五步骤 1. 对数据进行描述性分析，确定初始模型; 2. 估计模型中的未知参数; 3. 评估模型的拟合情况，并对初始模型进行必要的调整; 4. 检测模型的假设 (回归诊断); 5. 使用拟合的回归方程进行描述与预测 数据操作 数据类型/结构判断与转换 数据类型/结构判断与转换 • 回顾一下R有哪些数据类型/结构? • 给你一个数据，如何判断其数据类型/结构? • 如何把一个数据转化成所需要的数据类型/结构? 数据操作 补充 数据操作:排序相关 对于某一个序列，如 1.1，5.5, 3.3, 2.2, 4.4 • 排序 •秩 • 每一个是第几小(大)? • 哪个最大? • 哪个最小? 数据操作:补充 • 矩阵的行、列数、维数 • 矩阵的行列名称 • 子矩阵(分块矩阵) • 上三角阵，下三角阵 • 矩阵的加法、减法、乘法(常用三种乘法) • 谱分解/特征值分解(特征根，特征向量) • 奇异值分解(奇异值，左奇异向量，右奇异向量) • 主对角元、对角阵 数据操作 R中常见的数据结构与操作 列表 列表 g = "My First List" h = c(25, 26, 18, 39) j = matrix(1:10, nrow = 5) k = c("one", "two", "three") mylist = list(title = g, ages = h, j, k) mylist2 = c(title = g, ages = h, j, k) # output the 2nd element of mylist mylist[[2]] # output the element with label "ages" mylist[["ages"]] mylist$ages names(mylist) • 列表因为含有多种类型的变量(列表的元素甚至也可以是列表)，比较复杂， 使用的时候要小心。可以用str函数先查看一下列表的结构。 列表 • 假设下表是一个新冠肺炎病人的数据集。 • 我们从行的角度来看:每个病人的信息有多个，而且很多都不是一个类 型的。很多时候为了便于处理，我们希望把同一个病人的数据放在一起。 病人编号 入院时间 年龄 (PatientID) (AdmDate)(Age) 用药类型 病情 (Treatment) (Status) 1 01/22/2020 25 Remdesivir 危重 2 02/02/2020 34 Placebo 重 3 02/21/2020 28 Remdesivir 轻 4 01/08/2020 52 Remdesivir 危重 PatientID1 = 1 AdmDate1 = "01/22/2020" AdmDate1 = as.Date(AdmDate1, "%m/%d/%Y") Age1 = 25 Treatment1 = "Remdesivir" Status1 = "危重" 列表(将不同类型数据整合为一个对象) PatientData1a = c(PatientID1, AdmDate1, Age1, Treatment1, Status1) PatientData1a str(PatientData1a) PatientData1a[1] PatientData1b = list(PatientID1, AdmDate1, Age1, Treatment1, Status1) PatientData1b str(PatientData1b) PatientData1b[1] 列表(将不同类型数据整合为一个对象) PatientData1b = list(PatientID1, AdmDate1, Age1, Treatment1, Status1) PatientData1b str(PatientData1b) PatientData1b[1] PatientData1c = list(PatientID = PatientID1, AdmDate = AdmDate1, Age = Age1, Treatment = Treatment1, Status = Status1) PatientData1c str(PatientData1c) PatientData1c[1] 列表(如何访问列表内的数据) PatientData1c[1] PatientData1c["PatientID"] PatientData1c[[1]] PatientData1c[["PatientID"]] PatientData1c$PatientID 数据操作 R中常见的数据结构与操作 数据框 数据框(从列的角度来构造) 病人编号 入院时间 年龄 (PatientID) (AdmDate)(Age) 用药类型 病情 (Treatment) (Status) 1 01/22/2020 25 Remdesivir 危重 2 02/02/2020 34 Placebo 重 3 02/21/2020 28 Remdesivir 轻 4 01/08/2020 52 Remdesivir 危重 PatientID = c(1, 2, 3, 4) AdmDate = c("01/22/2020", "02/02/2020", "02/21/2020", "01/08/2020") AdmDate = as.Date(AdmDate, "%m/%d/%Y") Age = c(25, 34, 28, 52) Treatment = c("Remdesivir", "Placebo", "Remdesivir", "Remdesivir") Status = c("危重", "重", "轻", "危重") PatientData = data.frame(PatientID, AdmDate, Age, Treatment, Status) 数据框(从列的角度来构造) 病人编号 入院时间 年龄 (PatientID) (AdmDate)(Age) 用药类型 病情 (Treatment) (Status) 1 01/22/2020 25 Remdesivir Poor 2 02/02/2020 34 Placebo Improved 数据框( 3 02/21/2020 28 4 01/08/2020 52 Remdesivir Excellent Remdesivir Poor data frame):可以看成是一张表，其中每一列是同一种数据类型，而每 一行可以是不同数据类型。 PatientID = c(1, 2, 3, 4) AdmDate = c("01/22/2020", "02/02/2020", "02/21/2020", "01/08/2020") AdmDate = as.Date(AdmDate, "%m/%d/%Y") Age = c (25, 34, 28, 52) Treatment = c("Remdesivir", "Placebo", "Remdesivir", "Remdesivir") Status = c("危重", "重", "轻", "危重") PatientData = data.frame(PatientID, AdmDate, Age, Treatment, Status) PatientData = data.frame(PatientID, AdmDate, Age, Treatment, Status, stringsAsFactors = F) 注意最后两种构造方式中 stringsAsFactors选项 数据框元素的提取访问 PatientID1 = PatientData[1, 1] AdmDate1 = PatientData[1, 2] Age1 = PatientData[1, 3] Treatment1 = PatientData[1, 4] Status1 = PatientData[1, 5] class(PatientID1); str(PatientID1) class(AdmDate1); str(AdmDate1) class(Age1); str(Age1) class(Treatment1); str(Treatment1) class(Status1); str(Status1) 数据框变量的提取访问 # extract the 1st and 2nd column PatientData[1:2] PatientData[, 1:2] # extract the columns of "Treatment" and "Status" PatientData[c("Treatment", "Status")] # extract the variable Age PatientData$Age 数据框 几个问题: 1. 变量类型有些问题 str(PatientData) 用药类型是定性变量，而病情则是一个定序变量，这个信息没有记录。 病人编号 入院时间 年龄 (PatientID) (AdmDate)(Age) 用药类型 病情 (Treatment) (Status) 1 01/22/2020 25 Remdesivir Poor 2 02/02/2020 34 Placebo Improved 3 02/21/2020 28 Remdesivir Excellent 4 01/08/2020 52 Remdesivir Poor 数据框 • 用药类型是定性变量，而病情则是一个定序变量，这个信息没有记录。 PatientID = c(1, 2, 3, 4) AdmDate = c("01/22/2020", "02/02/2020", "02/21/2020", "01/08/2020") AdmDate = as.Date(AdmDate, "%m/%d/%Y") Age = c(25, 34, 28, 52) Treatment = c("Remdesivir", "Placebo", "Remdesivir", "Remdesivir") Treatment = factor(Treatment) Status = c("危重", "重", "轻", "危重") Status = factor(Status, ordered = T, levels = c("轻", "重", "危重")) str(PatientData1) PatientData1 = data.frame(PatientID, AdmDate, Age, Treatment, Status, stringsAsFactors = F) 注:一般来说，factor类型比字符类型节约空间，但是多了之后人很难记住每个后 面是什么意思。对于字符型向量，因子的水平默认依字母顺序创建。 数据框 • 可以在现有的数据基础上直接编码 PatientData2 = PatientData PatientData2$Treatment = factor(PatientData2$Treatment) PatientData2$Status = factor(PatientData2$Status, ordered = T, PatientData2 str(PatientData2) summary(PatientData2) levels = c("轻", "重", "危重")) 病人编号 入院时间 年龄 (PatientID) (AdmDate)(Age) 用药类型 病情 (Treatment) (Status) 1 01/22/2020 25 Remdesivir 危重 2 02/02/2020 34 Placebo 重 3 02/21/2020 28 Remdesivir 轻 4 01/08/2020 52 Remdesivir 危重 数据框 • 如果觉得每次都打"PatientData"太麻烦，可以用attach, detach, with函数。 • 以mtcars数据集为例，以下几种做法等价 summary ( mtcars$mpg ) () () plot plot mtcars$mpg, mtcars$disp mtcars$mpg, mtcars$wt attach ( summary mtcars ) ) ) ) ) ( mpg plot ( ( mpg, disp plot mpg, wt detach ( mtcars with { ) () () } ( mtcars, summary( plot plot mpg, disp mpg, wt mpg ) 数据清理 为什么要清理数据? 数据的清理 • 为什么要清理数据? • 原始数据可能存在错误、缺失等问题 • 实际工作中输入的数据未必以我们希望的方式存储与展示 • 对应于不同的分析工作，对于输入的数据有着不同的格式要求 • 综上，我们需要学会怎样把数据转换成符合我们要求的格式 数据清理 单数据集清理 单数据集清理 1. 数据行排序 2. 数据列换序 3. 改变元素的值 4. 处理缺失值 5. 删除含有缺失值的行(记录) 6. 修改现有的变量名 7. 对现有的变量类型进行重新设定 8. 选取部分数据(删除某些行列，保留某些行列) 9. 基于现有的数据创建新的数据(列) 单数据集清理 manager = c(1, 2, 3, 4, 5) date = c ("10/24/08", "10/28/08", "10/1/08", "10/12/08", "5/1/09") gender = c("M", "F", "F", "M", "F") age = c( 32, 45, 25, 39, 99) q1 = c(5 , 3, 3, 3, 2) q2 = c(4 , 5, 5, 3, 2) q3 = c(5 , 2, 5, 4, 1) q4 = c(5 , 5, 5, NA, 2) q5 = c(5 , 5, 2, NA, 1) leadership = data.frame(manager, date, gender, age, q1, q2, q3, q4, q5, stringsAsFactors = FALSE) leadership str(leadership ) leadership0 = leadership 单数据集清理(数据行排序) • 以哪些变量排序? • 如何排序?(升序还是降序) • 多个变量参加排序时如何设置优先级? # 按年龄(升序)排序 newdata <- leadership[order(leadership$age),] # 按年龄(降序)排序 newdata <- leadership[order(-leadership$age),] # 先按性别(升序)，再按年龄(升序)排序 attach(leadership) newdata <- leadership[order(gender, age),] detach(leadership) # 先按性别(升序)，再按年龄(降序)排序 attach(leadership) newdata <-leadership[order(gender, -age),] detach(leadership) 单数据集清理(数据列换序) • 对哪些列换序? • 如何排序?(按照数值还是名称) • 多个变量参加排序时如何设置优先级? # 按指定列的顺序换序 newdata = leadership[seq(from = 9, to = 1, by = -1)] newdata # 按指定列的名称换序 newdata = leadership[c("manager", "gender", "date", "age", "q1", "q2", "q3", "q4", "q5")] newdata 单数据集清理(缺失值处理) # 检查是不是缺失值 is.na(leadership) # 设置某些值为缺失值 leadership[age == 99, "age"] <- NA leadership # 分析中排除缺失值 x <- c(1, 2, NA, 3) y1 <- x[1] + x[2] + x[3] + x[4] y1 y2 <- sum(x) y2 y3 <- sum(x, na.rm=TRUE) y3 # 删除含有缺失值的行 newdata <- na.omit(leadership) newdata 单数据集清理(变量重新命名) • 把第二列变量名改为“testDate” names(leadership)[2] = "testDate" leadership 单数据集清理(对变量类型重新设定) # 对日期值做处理 leadership$date = as.Date(leadership$date, "%m/%d/%y") leadership str(leadership) 单数据集清理(数据抽取子集) • 保留部分变量 # 保留第6~9列 newdata <- leadership[, c(6:9)] # 保留指定列名的列 myvars <- c("q1", "q2", "q3", "q4", "q5") newdata <-leadership[myvars] • 删除部分变量 # 删除指定位置的变量 newdata <- leadership[,-c(8,9)] # 删除指定列名的变量 myvars <- names(leadership) %in% c("q3", "q4") newdata <- leadership[!myvars] 单数据集清理(数据抽取子集) • 保留部分观测 # 保留指定位置的观测 newdata <- leadership[1:3,] # 保留满足某些条件的观测 newdata <- leadership[leadership$gender=="M" & leadership$age > 30,] attach(leadership)
newdata <- leadership[gender=='M' & age > 30,]
detach(leadership)

单数据集清理(数据抽取子集)
• 保留部分数据(同时对行与列做筛选)
# 选择年龄不小于
35岁或小于24岁的观测，q1~q4 四个列
newdata
<- subset(leadership, age >= 35 | age < 24, select=c(q1, q2, q3, q4)) # 选择男性，年龄大于 25岁的观测，从gender到q4所有的列 newdata <- subset(leadership, gender=="M" & age > 25, select=gender:q4)

单数据集清理(数据抽取子集)
• 随机保留部分观测(抽样)
# 随机选择3个观测，无放回抽样
idx = sample(1:nrow(leadership), 3, replace = FALSE) newdata = leadership[idx, ]
newdata

单数据集清理(改变元素的值)
#创建agecat变量
leadership$agecat[leadership$age > 75] <- "Elder" leadership$agecat[leadership$age >= 55 & leadership$age <= 75] <- "Middle Aged" leadership$agecat[leadership$age < 55] <- "Young" # 更加紧凑的等价代码 leadership <- within(leadership,{ agecat <- NA agecat[age > 75] <- "Elder" agecat[age >= 55 & age <= 75] <- "Middle Aged" agecat[age < 55] <- "Young" }) 数据清理 多数据集清理 多数据集整理(数据集合并) • 列合并(直接合并) # combine two data sets horizontally leadership.A = leadership[, 1:2] leadership.B = leadership[, 3:9] leadership3 = cbind(leadership.A, leadership.B) leadership4 = data.frame(leadership.A, leadership.B) • 列合并(依据某一个变量) # merge two datasets by variable "manager" leadership.A = leadership[seq(from = 5, to = 1), 1:5] leadership.B = leadership[1:5, c(1, 6:9)] leadership5 = merge(leadership.A, leadership.B, by = "manager") 多数据集整理(数据集合并) • 行合并 leadership.A = leadership[1:3,] leadership.B = leadership[4:5,] leadership2 = rbind(leadership.A, leadership.B) • • 注意:此时两个数据框必须有相同的变量，但是顺序可以不同。如果 leadership.A 中有leadership.B没有的变量，会报错。 注意:如果两个都是 data frame，但是列名不同，也会报错。 数据清理 案例 案例:综合数据处理 我们要对以下学生评奖学金: 一共有数学、科学和英语三门考试。为了给所有学生确定一个单一的成绩衡量指标， 需要将这些科目的成绩组合起来。 你还想将前20%的学生评定为A，接下来20%的学生评定为B，依次类推。最后，你希望 按字母顺序对学生排序。 案例:综合数据处理 首先我们要输入数据 Student <- c("John Davis", "Angela Williams", "Bullwinkle Moose", "David Jones", "Janice Markhammer", "Cheryl Cushing", "Reuven Ytzrhak", "Greg Knox", "Joel England", "Mary Rayburn") Math <- c(502, 600, 412, 358, 495, 512, 410, 625, 573, 522) Science <- c(95, 99, 80, 82, 75, 85, 80, 95, 89, 86) English <- c(25, 22, 18, 15, 20, 28, 15, 30, 27, 18) roster <- data.frame(Student, Math, Science, English, stringsAsFactors=FALSE) 案例:综合数据处理 然后我们要处理成绩数据 # 标准化成绩 z <- scale(roster[,2:4]) # 求每个人的平均分 score <- apply(z, 1, mean) # 将平均分与原有数据合并 roster <- cbind(roster, score) # 计算分位数 y <- quantile(score, c(.8,.6,.4,.2)) # 给定分数等级 roster$grade[score >= y[1]] <- "A" roster$grade[score < y[1] & score >= y[2]] <- "B" roster$grade[score < y[2] & score >= y[3]] <- "C" roster$grade[score < y[3] & score >= y[4]] <- "D" roster$grade[score < y[4]] <- "F" 案例:综合数据处理 然后我们要处理姓名 # 将姓名按照空格分开 name <- strsplit((roster$Student), " ") # 提取姓 lastname <- sapply(name, "[", 2) # 提取名 firstname <- sapply(name, "[", 1) # 用分开的姓与名替代原来的姓名列 roster <- cbind(firstname,lastname, roster[,-1]) # 先按照姓(升序)，再按照名(升序)排序 roster <- roster[order(lastname,firstname),] # 输出结果 roster 数据清理 SQL语言 数据整理(SQL语言) • 数据库系统最流行的语言 • SQL = Structured Query Language • 前面学的语句仅可以在R中使用，在其他的系统中无法使用。 • 标准化的数据库操作语言:用于存取数据以及查询、更新和管理关系数 据库系统 • 广泛使用在各种数据库操作系统中，建议经常与数据打交道的同学花功 夫学习。 • 关系数据库(表与关系) • 收集所需要的数据 • 数据之间的关系清楚 • 尽量减少数据的冗余 • 常见数据库:Oracle，SQL-Server，SyBase，MySQL，...... 数据整理(SQL语言) • 数据查询语言(DQL, Data Query Language):也称为“数据检索语 句”，用以从表中获得数据，确定数据怎样在应用程序给出。 • 数据操作语言(DML, Data Manipulation Language):也称为动作查 询语言。用于添加，修改和删除表中的行。 • 事务处理语言(TPL, Transaction Process Language):它的语句能 确保被DML语句影响的表的所有行及时得以更新。 • 数据控制语言(DCL, Data Control Language):确定单个用户和用户 组对数据库对象的访问。 • 数据定义语言(DDL, Data Definition Language):在数据库中创建 新表或删除表;为表加入索引等。DDL包括许多与人数据库目录中获得 数据有关的保留字。它也是动作查询的一部分。 • 指针控制语言(CCL):用于对一个或多个表单独行的操作。 数据整理(SQL语言) 利用SQL检索数据 SQL检索语言的语法 Name Family Length beluga whale 15 dwarf shark 0.5 sperm whale 60 basking shark 30 humpback 50 whale shark 40 gray whale 50 blue whale 100 killer whale 30 mako shark 12 whale shark 40 SQL检索语言的语法 SELECT 表1.列1, 表1.列2, ... FROM Table1

;
• 从哪张表里查询哪些数据
• 这些数据按照什么格式展示(预处理)

查询marine中的数据(单张表)
* Inquire All Observations;
SELECT * FROM marine;
* Inquire Part Observations;
SELECT Name, Family
FROM marine;
• *表示所有列
• 变量之间要用逗号隔开
Name
Family
Length
beluga
whale
15
dwarf
shark
0.5
sperm
whale
60
basking
shark
30
humpback
50
whale
shark
40
gray
whale
50
blue
whale
100
killer
whale
30
mako
shark
12
whale
shark
40

查询marine中的数据(单张表)
* Inquire with Conditions;
SELECT *
FROM marine
WHERE Family != ”;
* Inquire with Conditions and Sorts;
SELECT *
FROM marine
WHERE Family != ”
ORDER BY Family, Length DESC;
• SQL语言中升序降序选项放在变量之后
Name
Family
Length
beluga
whale
15
dwarf
shark
0.5
sperm
whale
60
basking
shark
30
humpback
50
whale
shark
40
gray
whale
50
blue
whale
100
killer
whale
30
mako
shark
12
whale
shark
40

利用SQL合并数据

合并数据
class1
• 按列并(横向)，按行并(纵向)
• 注意:表class1中Alice不在表class2中，表class2中的Maggie不在表 class中
class2

合并数据(横向，内连接)
SELECT 表1.列1, 表1.列2, …, 表2.列1, 表2.列2, … FROM 表1, 表2
WHERE …
<其他语句>;
• 算法:
• 先构造笛卡尔积
• 再根据条件筛选处理
• 表之间要用逗号隔开

合并数据(横向，内连接)
SELECT c1.name, c1.sex, c1.age, c1.height, c2.weight FROM class1 AS c1, class2 AS c2
WHERE c1.name = c2.student_name;
• 注意:输出数据集中只有两个数据集共有的那些观测(Carol和James) class1 class2

合并数据(横向，外连接)
SELECT 表1.列1, 表1.列2, …, 表2.列1, 表2.列2, … FROM 表1
LEFT JOIN|RIGHT JOIN|FULL JOIN
表2
ON <连接条件> <其他语句>;

合并数据(横向，左连接)
SELECT *
FROM class1
LEFT JOIN class2
ON class1.name = class2.student_name;
class1
class2

合并数据(横向，右连接)
SELECT *
FROM class1
RIGHT JOIN class2
ON class1.name = class2.student_name;
class1
class2
目前sqldf包不支持!

合并数据(横向，全连接)
SELECT *
FROM class1
FULL JOIN class2
ON class1.name = class2.student_name;
class1
class2
目前sqldf包不支持!

合并数据(纵向)
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d

合并数据(纵向)
SELECT * FROM 表A <其他语句>
连接方式 SELECT * FROM 表B <其他语句>;
• 连接方式:
• EXCEPT:表A中去掉重复的行，在表A但不在表B中的行。 • INTERSECT:选择在表A中不重复的，同时在表B中的行 • UNION:在表A中或在表B中，且不重复的所有行
• OUTER UNION:同时查看表A和表B中的所有列
• ALL:保留表A中重复的行
• CORR:根据列的名称进行合并(否则只看列类型是否相同)，删除所有不是同时
在两个表中的列。

合并数据(纵向，EXCEPT)
SELECT * FROM A EXCEPT
SELECT * FROM B;
取𝐴\B,重复的观测不保留
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d

合并数据(纵向，EXCEPT ALL)
SELECT * FROM A EXCEPT ALL SELECT * FROM B;
取𝐴\B,重复的观测保留 目前sqldf包不支持!
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d

合并数据(纵向，EXCEPT CORR)
SELECT * FROM A EXCEPT CORR SELECT * FROM B;
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d
取A与B同名的列，再求𝐴\B,重复的观测不保留 目前sqldf包不支持!

合并数据(纵向，EXCEPT ALL CORR)
SELECT * FROM A EXCEPT ALL CORR SELECT * FROM B;
取A与B同名的列，再求𝐴\B,重复的观测保留 目前sqldf包不支持!
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d

合并数据(纵向，INTERSECT)
SELECT * FROM A INTERSECT
SELECT * FROM B;
求𝐴 ∩ B,重复的观测不保留
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d

合并数据(纵向，INTERSECT ALL)
SELECT * FROM A INTERSECT ALL SELECT * FROM B;
求𝐴 ∩ B,重复的观测保留 目前sqldf包不支持!
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d

合并数据(纵向，INTERSECT CORR)
SELECT * FROM A INTERSECT CORR SELECT * FROM B;
取A与B同名的列，再求𝐴 ∩ B,重复的观测不保留 目前sqldf包不支持!
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d

合并数据(纵向，INTERSECT ALL CORR)
SELECT * FROM A INTERSECT ALL CORR SELECT * FROM B;
取A与B同名的列，再求𝐴 ∩ B,重复的观测保留 目前sqldf包不支持!
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d

合并数据(纵向，UNION)
SELECT * FROM A UNION
SELECT * FROM B;
求𝐴 ∪ B,重复的观测不保留
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d

合并数据(纵向，UNION ALL)
SELECT * FROM A UNION ALL
SELECT * FROM B;
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d
取A与B同名的列，再求𝐴 ∪ B,重复的观测保留

合并数据(纵向，UNION CORR)
SELECT * FROM A UNION CORR SELECT * FROM B;
求𝐴 ∪ B,重复的观测不保留 目前sqldf包不支持!
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d

合并数据(纵向，UNION ALL CORR)
SELECT * FROM A UNION ALL CORR SELECT * FROM B;
取A与B同名的列，再求𝐴 ∪ B,重复的观测保留 目前sqldf包不支持!
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d

合并数据(纵向，OUTER UNION)
SELECT * FROM A OUTER UNION SELECT * FROM B;
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d
目前sqldf包不支持!

合并数据(纵向，OUTER UNION CORR)
SELECT * FROM A OUTER UNION CORR SELECT * FROM B;
表A 表B
ID
X
1
a
1
b
1
c
1
a
2
b
ID
Y
1
b
2
b
3
c
3
c
4
d
目前sqldf包不支持!

R作图

R中绘图书籍

R中绘图函数的种类
在R中有两种绘图函数 • 高级绘图函数
• 创建一个新的图形
• 低级绘图函数
• 在现有的图形上添加元素,对原有图进行补充

高级绘图函数
函数名 功能
plot(x) 以x的元素值为纵坐标、以序号为横坐标绘图 plot(x,y) x与y的二元散点图
pie(x) 饼图
boxplot(x) 盒形图(也称箱线图)
hist(x) x的频率直方图
barplot(x) x的值的条形图
qqnorm(x) 正态分位数-分位数图

高级绘图函数的选项
选项
add=FALSE axes=TRUE type=”p”
xlim=,ylim=
xlab=,ylab=
main=
sub=
功能 如果是TRUE，叠加图形到前一个图上(如果有的话) 如果是FALSE，不绘制轴与边框
指定图形的类型，”p”:点，”l”:线，”b”:点连线，”o”: 同上，但是线在点上，”h”:垂直线，”s”:阶梯式，垂直 线顶端显示数据，”S”:同上，但是垂直线底端显示数据
指定轴的显示范围
坐标轴的标签
主标题
副标题

低级绘图函数
函数名
points(x,y) lines(X,Y) text(x,y,labels,…) abline(a,b) abline(h=y) abline(v=x) abline(lm.obj) legend(x,y,legend)
title()
axis(side,vect)
功能
添加点
添加线
加标记，在(x,y)处添加用labels指定的文字
加直线，绘制斜率为b和截距为a的直线
加直线，在纵坐标y处画水平线
加直线，在横坐标x处画垂直线
加直线，画出lm.obj确定的回归线
加注释，在点(x,y)处，说明内容由legend给 定
加标题，也可添加一个副标题
加坐标轴

图形的元素与对应的参数
main
ylim=c(110,170)
type=”b”
ylab
axes
xlim=c(50,80)
xlab sub

图形元素示例(关于边框)
plot(women
$height,women$weight)
plot(women
$height,women$weight,type=”p”,
main=”Women Data”,sub=”Ten Women 30-39″,
,font.lab=1,xlab=”height”,ylab=”weight”,
xlim=c(50,80),ylim=c(110,170))
plot(women
$height,women$weight,type=”p”,
main=”Women Data”,sub=”Ten Women 30-39″,
,xlab=”height”,ylab=”weight”,
font.lab=1,axes=F)

图形元素示例(关于线的类型)
par(mfrow=c(2,3)) plot(women$height,women$weight,xlab=””,ylab=””,font.lab=1,
type=”p”,main=”type is p”) plot(women$height,women$weight,xlab=””,ylab=””,font.lab=1,
type=”l”,main=”type is l”) plot(women$height,women$weight,xlab=””,ylab=””,font.lab=1,
type=”b”,main=”type is b”) plot(women$height,women$weight,xlab=””,ylab=””,font.lab=1,
type=”o”,main=”type is o”) plot(women$height,women$weight,xlab=””,ylab=””,font.lab=1,
type=”h”,main=”type is h”) plot(women$height,women$weight,xlab=””,ylab=””,font.lab=1,
type=”s”,main=”type is s”) par(mfrow=c(1,1))
dev.off()

图形元素示例(关于线的类型)
“p”:点
“l”:线
“b”:点连线 “o”:同上，但是线在点上 “h”:垂直线 “s”:阶梯式，

其他常用图形元素
选项
adj
cex
col
font lwd
lty
pch
pty
xaxt, yaxt
功能
控制关于文字对齐方式(0-左对齐，0.5-居中对齐，1-右 对齐)
符号和文字大小,cex.axis,cex.lab,cex.main,cex.sub 颜色,col.axis,col.lab,col.main,col.sub
文字字体(1.正常，2.粗体，3.斜体，4.粗斜体)， font.axis, font.lab,font.main,font.sub
线的宽度
连线的线型(1.实线，2.虚线，3.点线，4.点虚线，5。长 虚线，6.双虚线)
绘图符号的类型(1到25的整数) 绘图区域类型
如果xaxt=”n”,设置x轴不显示,如果yaxt=”n”,设置y轴不 显示

其他常用图形元素(边界)

其他常用图形元素(边界)
选项
mar
mai mgp
mfcol, mfrow
功能 图形边空的大小c(bottom,left,top,right),缺
省值为c(5.1,4.1,4.1,2.1)，单位是英寸 同上，单位是文本行数
轴标题，轴标签和轴线到坐标轴的距离, 单位是文本行
c(nr,nc)的向量，前者按列分隔绘图窗口, 后者按行分隔绘图窗口
par(): 设置全局绘图参数

其他常用图形元素(大小等)

其他常用图形元素(大小等)
par(mgp = c(2, 1, 0.2), mar = c(4, 3, 2, 1)) ###修改缺省值
X=1:3; Y=1:3
plot(X, Y, col = “blue”,
pch = 16, #字体大小
cex = c(0.5, 1, 2), #绘图符号大小
cex.axis = 1.1, #坐标轴刻度数字大小
cex.lab = 1.2, #坐标轴标签大小
cex.main = 1.4, #标题文字大小
cex.sub = 1.3, #副标题文字大小
font.axis = 3, #坐标轴刻度数字字体，斜体
font.lab = 2, #坐标轴标签字体，粗体
font.main = 4, #标题字体，粗斜体
)
font.sub = 1, #副标题字体，正常
main = “Example”,
xlim = c(1, 12), ylim = c(1, 12)
text(X, Y, adj=0, labels=paste(“cex=”, c(0.5,1,2)), col=”red”, cex=1.2, font=3)

其他常用图形元素(符号)
par(mgp = c(1.6, 0.6, 0), mar = c(3, 3, 2, 1))
pch_type = 1:25
X = 1:25;Y = rep(6, 25)
plot(X, Y, col = 1, pch = pch_type, cex = 2,
main = “pch:1-25”, font.lab = 2)
text(X, Y, adj = -0.5,labels=paste(“pch=”,pch_type),srt=-90)

其他常用图形元素(颜色)
# 1, 颜色名称:white, black, red,…
colors()[1:10]
# 2″数字代号:1, 2, 3, 4, 5, 6, 7, 8,…
pch_type=1:25
X=1:25;Y=rep(6,25)
plot(X,Y,col=1:25,pch=pch_type,cex=2,main=”pch:1-25″,font.lab=2)
text(X,Y ,col=1:25,adj=-0.5,labels=paste(“pch=”,pch_type),srt=-90)

其他常用图形元素(颜色)
SoftDrink = read.csv(“SoftDrink.txt”,header=F)
SoftDrink = SoftDrink[SoftDrink!=””]
# check how many types of soft drink
unique(SoftDrink)
table(SoftDrink)
prop=prop.table(table(SoftDrink))
barplot(prop,col=rainbow(5))
slices.1 <- as.numeric(table(SoftDrink)) labels.1 <- names(table(SoftDrink)) pie(slices.1, col=rainbow(5), labels = paste(labels.1,":",prop,"%"), main="Pie Chart of Soft Drinks", cex.lab=0.5) 其他常用图形元素(线型) #lty 控制连线的类型: #1: 实线，2: 虚线，3: 点线， #4: 点虚线，5: 长虚线，6: 双虚线 #lwd 控制线的宽度: #1为正常大小，<1 为减小宽度，>1为增加宽度
par(mgp=c(1.6,0.6,0),mar=c(3,3,2,1))
X=1:10;Y=1:10
line_type=c(1:6)
line_width=seq(from=0.6,by=0.6,length=6)
plot(X,Y,col=”red”,pch=16,type=”n”,
main=”line type and line width”,
font.lab=2,cex.lab=1.2)
abline(h=3:8,lty=line_type,col=1:6,lwd=line_width)
text(3:8,2.8:7.8,adj=0,
labels=paste(“lty=”,line_type,”,lwd=”,line_width),
col=”blue”,cex=1.2,font=3)

小综合(点线混合)
plot(women[,1], women[,2], pch=16, col=”red”,
xlab=”weight”,ylab=”height”, type=”p”,font.axis=2,
font.lab=2,cex.lab=1.5)
lm(women[,2]~women[,1])
x=seq(58,72,by=0.1)
y=-87.52+3.45*x
lines(x,y,col=”green”,lwd=2)
a=women[,1]; b=-87.52+3.45*a
points(a,b,pch=15,col=”blue”)
legend(“topleft”,pch=c(16,-1,15),lty=c(-1,1,-1),
col=c(“red”, “green”,”blue”),
legend=c(“original points”,”lines”,”fitted points”))

其他常用图形元素
函数名 功能
pairs(m) 如果x是矩阵或是数据框，作x的各列之间的二元图 coplot(x~y|z) 关于z的每个数值(或数值区间)绘制x与y的二元图 matplot(x,y) 二元图，其中x的第一列对应y的第一列，依次类推 dotchart(x) 点图
image(x,y,z) 三维图形 contour(x,y,z) 等值线
persp(x,y,z) 三维图形的表面曲线 draw.circle() 韦恩图
pie3D 三维饼图
其他图形包:gplots,ggplot2,maps,…

统计软件
第八讲:R语言初步5
黄达 复旦大学管理学院统计学系 dahuang@fudan.edu.cn

编程 控制流

控制流:基本概念
• 语句(statement)是一条单独的R语句或一组复合语句(包含在花括 号{ } 中的一组R语句，使用分号分隔);
• 条件(cond)是一条最终被解析为真(TRUE)或假(FALSE)的表达 式;
• 表达式(expr)是一条数值或字符串的求值语句;
• 序列(seq)是一个数值或字符串序列。
•
•
• /选择 •
三种控制流
顺序
条件
循环

控制流:条件执行
• 结构
语法为:
if-else
控制结构if-else在某个给定条件为真时执行语句。也可以同时在条件为假时执
行另外的语句。
if(cond) statement
if(cond) statement1 else statement2
# 如果grade是字符型，转成factor型
if(is.character(grade)) grade <- as.factor(grade) # 如果grade不是factor型，转成factor型 # 否则输出"Grade already is a factor." if(!is.factor(grade)){ grade <- as.factor(grade) }else{ print("Grade already is a factor.") } 控制流:条件执行 • 结构 ifelse ifelse 结构是if-else结构比较紧凑的向量化版本。 ifelse (cond, statement1, statement2) 语法为: # 如果成绩大于 0.5，显示"Passed"，否则显示"Failed" ifelse (score>0.5,print(“Passed”),print(“Failed”))
outcome
<- ifelse(score>0.5, print(“Passed”),print(“Failed”))

控制流:条件执行
• 结构
switch
switch根据一个表达式的值选择语句执行。
语法为:
switch(expr, …)
feelings <- c("sad", "afraid") for (i in feelings){ print( switch(i, happy = "I am glad you are happy", afraid = "There is nothing to fear", sad = "Cheer up", angry = "Calm down now" ) ) } 控制流:重复与循环 • for结构:重复地执行一个/组语句，直到某个变量的值不再包含在序列seq中为止。 语法为: # 打印1~10 for(i in 1:10) for (var in seq) statement print(i) # 打印"Round 1"~"Round 10" for(i in 1:10){ string = paste("Round",i) print(string) } 控制流:重复与循环 • while结构:重复地执行一个/组语句，直到条件不为真为止。语法为: # 打印"Round 10"~"Round 1" i = 10 while(i > 0){
print(paste(“Round”,i))
i=i–1
}
while
(cond) statement

控制流:重复与循环
• repreat结构:重复地执行一个/组语句，直到条件不为真为止。语法为:
repeat{
statements…
if(condition){
}
break
}
# 打印”Round 10″~”Round 1″
i = 10
repeat{
print(paste(“Round”,i))
i <- i – 1 if(i <= 0){ } break } 控制流:重复与循环 • break语句 程序从循环中跳出。 • next语句 程序忽略其余语句，直接进入下一次循环。 • 备注: 在处理大数据集中的行和列时，R中的循环可能比较低效费时。 只要可能，最好联用 R中的内建数值/字符处理函数和apply族函数。 补充:节约循环时间 补充:节约循环时间 • 由于R是解释型语言，所以在执行循环语句时，每次都要把代码转成机 器执行指令，因此效率低下。 • 生成从1加到N，保留每一个结果，并计算时间 • 比较一下下列三个程序: • 尽量多用内置函数或者是向量计算，避免循环 N = 50000 N = 50000 N = 50000 time0 = Sys.time() x = rep(0,N) for(i in 1:N){ x[i] = sum(1:i) } time1 = Sys.time() print(time1-time0) time0 = Sys.time() time0 = Sys.time() x=1 x = 1:N for(i in 2:N){ x = cumsum(x) x = c(x,sum(1:i)) time1 = Sys.time() } print(time1-time0) time1 = Sys.time() print(time1-time0) 编程 自定义函数 编程:自定义函数 • 很多时候R自己带的函数不能满足我们的需求，这个时候需要自己编写， 构造一个函数。 • R中自定义函数的一般形式是: myfunction <- function(arg1, arg2, ...){ statements return(object) } • 函数中的对象只在函数内部使用。 • 返回对象的数据类型是任意的，从标量到列表皆可。(一次可以返回 很多信息!) 编程:自定义函数 • 假设你想编写一个函数 • 输入:一组数值 • 任务: • 数据对象的集中趋势和散布情况。此函数可以选择性地给出参数统计量(均值 和标准差)和非参数统计量(中位数和绝对中位差)。 • 输出 • 结果以一个含名称列表的形式给出。 • 用户可以选择是否自动输出结果。除非另外指定，否则此函数的默认行为应当 是计算参数统计量并且不输出结果。 编程:自定义函数 mystats <- function(x, parametric=TRUE, print=FALSE) { if (parametric ){ center } else <- mean(x); spread <- sd(x) { center <- median(x); spread <- mad(x) } if (print & parametric) { cat ("Mean=", center, "\n", "SD=", spread, "\n") } else if (print & !parametric) { cat ("Median=", center, "\n", "MAD=", spread, "\n") } } result <- list(center=center, spread=spread) return (result) set.seed (1234) x <- rnorm (500) y <- mystats (x) y <- mystats (x, parametric=FALSE, print=TRUE) 编程:自定义函数 • 编写一个函数 • 输入:日期格式 • 功能:按照指定格式输出当前的日期 • 输出:日期(含指定的格式) • 备注:设置默认值(在函数声明中为参数指定值) mydate } <- function(type="long") { switch (type, long = format(Sys.time(), "%A %B %d %Y"), short = format(Sys.time(), "%m-%d-%y"), cat (type, "is not a recognized type\n")) mydate ("long") mydate mydate ("short") () mydate ("medium") 一个综合案例:图形的SVD分解 图形的SVD分解 图形的SVD分解 图形的SVD分解 rawimg = readJPEG("pic0.jpg"); lst = factorize(rawimg, 100); neig = c(1, 5, 20, 50, 100); for(i in neig){ m = recoverimg(lst, i); writeJPEG (m, sprintf("svd_%d.jpg", i), 0.95); } 本程序由邱怡轩博士提供，特此致谢! 图形的SVD分解 # By Yixuan Qiu (COS) library ("rARPACK") library ("jpeg") factorize = function(m, k){ r = svds(m[,,1], k) g = svds(m[,,2], k) b = svds(m[,,3], k) return (list(r = r, g = g, b = b)); } 本程序由邱怡轩博士提供，特此致谢! 图形的SVD分解 recoverimg = function(lst, k){ recover0 = function(fac, k){ dmat = diag(k); m = fac$u[, 1:k] %*% dmat %*% t(fac$v[, 1:k]); m[m < 0] = 0; m[m > 1] = 1;
return(m);
}
diag(dmat) = fac$d[1:k];
r=
recover0(lst$r, k);
g=
recover0(lst$g, k);
b=
recover0(lst$b, k);
m=
array(0, c(nrow(r), ncol(r), 3));
m[, ,
1] = r;
m[, ,
2] = g;
m[, ,
3] = b;
return
(m);
}
本程序由邱怡轩博士提供，特此致谢!

统计软件
第十讲:描述统计1
黄达 复旦大学管理学院统计学系 dahuang@fudan.edu.cn

数据中有哪些信息?
怎么从数据中提取信息?

数据中的信息
• 离散型数据
源自Statistics for Business and Economics (11th Edition)

数据中的信息
• 离散型数据
源自Statistics for Business and Economics (11th Edition)

数据中的信息
• 离散型数据
源自Statistics for Business and Economics (11th Edition)

数据中的信息
• 离散型数据
源自Statistics for Business and Economics (11th Edition)

数据中的信息
• 离散型数据
源自Statistics for Business and Economics (11th Edition)

数据中的信息
SoftDrink = read.csv(“SoftDrink.txt”,header=F,stringsAsFactors=F) SoftDrink = SoftDrink[SoftDrink!=””]
unique(SoftDrink) # 列出所有的软饮料(去重复的) table(SoftDrink)
prop.table(table(SoftDrink)) barplot(prop.table(table(SoftDrink)))
slices.1 <- as.numeric(table(SoftDrink)) labels.1 <- names(table(SoftDrink)) pie(slices.1, labels = labels.1, main="Pie Chart of Soft Drinks") 数据中的信息 pie(x, labels = names(x), edges = 200, radius = 0.8, clockwise = FALSE, ...) require (grDevices) pie(rep (1, 24), col = rainbow(24), radius = 0.9) pie.sales <- c(0.12, 0.3, 0.26, 0.16, 0.04, 0.12) names(pie.sales ) <- c("Blueberry", "Cherry", "Apple", "Boston Cream", "Other", "Vanilla Cream") pie(pie.sales ) # default colours pie(pie.sales , col = c("purple", "violetred1", "green3", "cornsilk", "cyan", "white")) pie(pie.sales , col = gray(seq(0.4, 1.0, length = 6))) pie(pie.sales , density = 10, angle = 15 + 10 * 1:6) pie(pie.sales , clockwise = TRUE, main = "pie(*, clockwise = TRUE)") 数据中的信息 barplot (height, width = 1, space = NULL, names.arg = NULL, legend.text = NULL, beside = FALSE, horiz = FALSE, density = NULL, angle = 45, col = NULL, border = par("fg"), ...) require (grDevices) # for colors tN <- table (Ni <- stats::rpois(100, lambda = 5)) r <- barplot (tN, col = rainbow(20)) #- type = "h" plotting *is* ' bar'plot lines(r, tN, type = "h", col = "red", lwd = 2) 数据中的信息(表格) • 连续型数据 源自Statistics for Business and Economics (11th Edition) 数据中的信息(表格) • 连续型数据 源自Statistics for Business and Economics (11th Edition) 数据中的信息(表格) • 连续型数据 源自Statistics for Business and Economics (11th Edition) 数据中的信息(图形) • 连续型数据 源自Statistics for Business and Economics (11th Edition) 数据中的信息(图形) • 连续型数据 v.s. 离散型数据 源自Statistics for Business and Economics (11th Edition) 数据中的信息(表格) • 连续型数据:累计频率表(回忆CDF) 源自Statistics for Business and Economics (11th Edition) 数据中的信息(图) • 连续型数据:累计频率图(回忆CDF) 源自Statistics for Business and Economics (11th Edition) 数据中的信息 Audit = read.table("Audit.txt",header=F) Audit = as.vector(as.matrix(Audit)) region = c(10,14,19,24,29,34) Table.Freq = table(cut(Audit,region)) Table.Freq Table.Freq2 = prop.table(Table.Freq) x.label = "Audit Time (days)" y.label = "Frequency" # plot histogram (with warning) hist(Audit,breaks=region,xlim=c(0,35),freq=TRUE,ylab=y.label,xlab=x.label) # plot histogram (without warning) hist(Audit,breaks=c(9,14,19,24,29,34),xlim=c(0,35),freq=TRUE,ylab=y.label,xlab=x.lab el) cdf.Audit = cumsum(Table.Freq2) y.label = "Cummulative Frequency" plot(cdf.Audit,type="b",ylab=y.label,xlab=x.label) 数据中的信息 数据中的信息 library("readxl") SHStock = read_excel("SH000001.xlsx") n = nrow(SHStock) SHA = unlist(SHStock[, 4]) rate.A = log(SHA[2:n]/SHA[1:(n-1)]) hist(rate.A, xlim = c(-0.4, 0.4), breaks = seq(-0.8,0.8,0.01), main = "上证指数收益率(基于开盘价)") SHB = unlist(SHStock[, 7]) rate.B = log(SHB[2:n]/SHB[1:(n-1)]) hist(rate.B, xlim = c(-0.4, 0.4), breaks = seq(-0.8,0.8,0.01), main = "上证指数收益率(基于收盘价)") 数据中的信息(抢红包) • 你应该先抢还是后抢? • 先抢比后抢会多抢到钱么?还是后抢比先抢会抢 到更多的钱? • 不同的人数抢红包，谁更有可能抢到“最佳手 气”? 源自《一席》毕啸天演讲《一本正经 胡说八道》 数据中的信息(抢红包) • 5人，50元红包，150次 源自《一席》毕啸天演讲《一本正经 胡说八道》 数据中的信息(抢红包) • 5人，50元红包，5000万次 源自《一席》毕啸天演讲《一本正经 胡说八道》 数据中的信息(抢红包) • N个人抢N个红包，抢到最佳手气的概率 源自《一席》毕啸天演讲《一本正经 胡说八道》 数据中的信息(图与表小结) 其他常见的图 数据中的信息(图) • 网络图 源自Statistical Inference 数据中的信息(图) • 网络图 源自网络 数据中的信息(图) • 词云图 源自网络 数据中的信息(图) • 词云图 源自网络 数据中的信息(图) • 泡泡图 • http://assets.dtcj.com/visualization/metro/metro_in_a_day.html?from=timel ine&isappinstalled=0 源自网络 数据中的信息(图) • 热力图 源自百度 数据中的信息(图) • 热力图 前Boston Red Sox棒球队传奇巨星 Ted Williams一次比赛的照片。 源自网络 数据中的信息(图) • 热力图 源自网络 数据中的信息(图) • 时间序列图 源自网络 数据中的信息(图) • 一张好的图应该包含的信息: • 标题 • 坐标含义 • 坐标的方向、单位 • 数据所表达的意思 • 数据来源 • 不能遮掩、误导信息 建议 拿到数据后先不要急着跑模型，先做一些 表和图，这样会对数据有第一手的认识。 画图的实际应用 实际应用:识别僵尸粉 • 幂律(厚尾)分布与查找“僵尸粉” 实际应用:发改委“打飞机” • 2009年3月25日，发改委决定将汽、柴油价格每吨分别提高290元和 180元;当天，美国空军一架正在执行测试飞行F-22“猛禽”战机 加州爱德华兹空军基地以北六英里的地方坠毁。 • 2009年6月1日，发改委发布调价通知，上涨油价;当天14时,一架 载有228人的法航空客A330起飞不久后与地面失去联系。 • 2009年6月30日，发改委再次发布调价通知，上涨油价;一架载有 154人的客机在从也门前往科摩罗的途中坠毁。 • 2009年7月15日，发改委就成品油价格问题发表说明，称价格未调 整到位。当日，伊朗里海航空公司的一架客机在该国西北部城市加 兹温附近的村庄坠毁。 • 2009年7月24日，近期部分城市提高水价引发社会广泛关注，发改 委当天表示，要合理把握水价调整的力度和时机。当天，一架由德 黑兰飞往马什哈德的伊朗阿利亚航空公司的客机在马什哈德国际机 场降落时滑出跑道。 源自 http://blog.sina.cn/dpool/blog/s/blog_6fea8c730101lql k.html?vt=4 实际应用:发改委“打飞机” • 2009年8月4日，国家发展和改革委员会能源研究所原所长周大 地4日表示，随着天然气供应量的逐渐增加，价格的提高这也 是不可避免的。当天，一架曼谷航空公司飞机于下午三点十六 分在泰国南部苏梅岛着陆时失事。 • 2009年9月3日，发改委官员释疑油价上调推迟:需考虑更多因 素。一架直升机在美国密西西比州首府杰克逊市附近坠毁，砸 毁地面一处无人居住的房子，但没造成地面人员伤亡。机上载 有美国联 邦航空管理局两名监察员，其中一人已死亡，另外 一人受伤。 • 2009年9月29日，汽、柴油价格每吨均下调190元。降的不多不 少，只伤飞机:当天，据德国《图片报》9月29日报道，近日， 一架载着80名乘客、 飞往科索沃首府普里什蒂纳的波音737- 700客机，刚刚从杜塞尔多夫机场起飞，就遭到一群欧掠鸟的 “袭击...... 源自 http://blog.sina.cn/dpool/blog/s/blog_6fea8c730101lql k.html?vt=4 实际应用:发改委“打飞机” • 2009 年11月9日，发改委再次宣布柴油每吨提价428。结果晚上肯 尼亚果然又掉飞机了，新华网报道，一架小型货机9日在肯尼亚首 都内罗毕一个机场坠毁。 • 2010年3月25日，发改委副主任彭森25日表示，2010年将重点推进 居民用水、用电阶梯式价格改革，年内还要推进资源税改革，研究 制定开征环境税的改革方案。一架大马皇家空军部队的飞机，25日 上午在北方大学体育大楼附近的山区坠毁。 • 2010年4月1日，发改委宣布国内航线旅客燃油附加费率的通知。美 国第五舰队发表声明称，一架载有4名机组人员的美海军E-2C“鹰 眼”预警机31日在阿拉伯海坠毁。 • 2010年4月2日，据武汉晨报4月2日报道，国家发改委相关司局近日 召开会议，酝酿重新启动电网输配分开改革。当地时间4月2号，巴 西南部一架小型飞机在进行飞行表演时意外坠地，并爆炸起火。 • 2010年4月13日，发改委宣布97#汽油今起涨至6.92元。香港机场一 客机着陆时轮胎引发火警;印尼机场一客机降落后滑出跑道坠河; 源自 http://blog.sina.cn/dpool/blog/s/blog_6fea8c730101lql k.html?vt=4 实际应用:发改委“打飞机” • 2010年4月16日，发改委暗示将上调资源价格。当日，印度自制火 箭发射失败。 • 2010年4月21日，发改委公布全国36个大中城市房价监测数据,数据 显示,3月份全国房价继续上涨。当日，哥伦比亚2直升机相撞坠毁。 • 2010年5月12日，发改委表示，预计5、6月份煤炭价格总体将保持 稳定，后期可能小幅上涨。同日，利比亚一架客机坠毁，机上105 人团灭。 • 2010年10月25日，25日晚9时，发改委宣布自2010年10月26日0时起 汽柴油价格每吨分别提高230元和220元。26日，位于南非首都比勒 陀利亚的温德博恩机场发生飞机相撞事故，造成1人重伤，2人轻伤。 • 2010年12月21日，发改委公布上调汽油柴油价格。美国国民警卫队 一架载有6人的直升飞机在波多黎各坠毁。 源自 http://blog.sina.cn/dpool/blog/s/blog_6fea8c730101lql k.html?vt=4 实际应用:发改委“打飞机” • 2011年4月6日，发改委决定，自2011年4月7日零时起，将汽柴 油销售价格每吨分别上调500元和400元，相当于90号汽油和0 号柴油每升分别提高0.37元和0.34元。4月6日俄空军一架苏- 27战斗机当天执行例行飞行任务时在滨海边疆区失事。同日， 美军一架F/A-18超级“大黄蜂”战机6日中午在加利福尼亚州 中部地区坠毁。另有土尔其一架直升机坠毁。 • 2012年11月16日 ，据国家发展改革委通知，决定自今日零时 起将汽、柴油价格每吨分别降低310元和300元，这是今年国内 成品油价第四次下调。美国WMBB网站报道，当地时间11月15日 15时，美国泰恩代尔空军基地的一架F-22A战斗机，在美国本 土佛罗里达州的98号公路附近坠毁。F-22A的飞行员在飞机坠 毁前，安全弹射，此次事故目前并未造成人员伤亡。 源自 http://blog.sina.cn/dpool/blog/s/blog_6fea8c730101lql k.html?vt=4 实际应用:发改委“打飞机” 我们庄重承诺: 绝不首 先使用 发改委 源自网络 实际应用:发改委“打飞机” • 在2005年3月23日至2011年4月7日期间(2206天)，发改委上调油价16次，期 间飞机失事126次。平均每17.5天有一架飞机失事。 图有异常，必有妖孽 数据中的信息(图) 源自《上市公司连续两年亏损就应该被ST吗》， 姜国华，王汉生，《经济研究》2005年第3期 数据中的信息(图) 源自网络 常见图表忽悠法 (己之不欲，勿施于人) 本部分来自:http://simplystatistics.org/2012/11/26/the-statisticians-at-fox-news-use-classic- and-novel-graphical-techniques-to-lead-with-data/ 忽悠法1 修改横轴长度 源自网络 忽悠法2 坐标变换 源自网络 忽悠法2:或许不是忽悠 • 时间序列图 源自百度 忽悠法3 坐标截断 源自网络 忽悠法3 坐标截断 源自网络 忽悠法4 选择部分数据 源自网络 忽悠法5 修改数字 源自网络 结论 • 光看图是不够的，容易被误导，还得看数。 描述统计 描述统计 • 如果有很多数，比如说10000个数:𝑥 , ⋯ , 𝑥 1 10000 • 显然，我们要是每个数都看的话，是不现实的。 • 所以，我们需要将这些数浓缩成若干个指标(描述统计量)，把这堆数 中的重要信息提取出来。 • 拿到一个数据集后我们想做什么? • 数据内含有信息，我们想获取数据中的信息。 • 要获取所有的信息，最好的办法是所有的数据都看一遍。数据量很大的时候我 们无法做到。 • 我们希望能够提取出一些少量但是重要的信息出来 • 讨论: • 哪些信息重要?(举例说明) • 哪些指标对应着这些信息? 预备知识 • 假定有比率数据 𝑥 ,𝑥 ,⋯,𝑥 12𝑛 •排序后的数据 𝑥(1),𝑥(2),⋯,𝑥(𝑛) •数据的和:𝑆 =σ𝑛 𝑥 =𝑥 +𝑥 +⋯+𝑥 𝑛𝑖=1𝑖12 𝑛 • 有时候简写成σ𝑥𝑖 单变量数据的描述统计 如何刻画“集中位置” 描述统计(众数) • 平均涉及到四则运算 • 定性数据无法进行四则运算，不存在平均数 • 问:如何解决? • 答:众数?(出现次数最多的那类) • 例子:重大事项的投票 • 2/3以上出席，出席者1/2以上通过 描述统计(简单算术平均) • 假定现在你手里有𝑘支股票 • 每只股票你只有1股 • 每只股票1元/股 • 股票收益率依次为:𝑥 ,𝑥 ,⋯,𝑥 • 问题:你如何衡量你的收益? 𝑥ҧ=𝑥1+𝑥2+⋯+𝑥𝑘 =1σ𝑘 𝑥 𝑘 𝑘 𝑖=1 𝑖 12𝑘 描述统计(加权算术平均) • 假定现在你手里有𝑘支股票 • 第𝑖只股票，有𝑛𝑖股 • 每只股票1元/股 • 股票收益率依次为:𝑥 ,𝑥 ,⋯,𝑥 12𝑘 • 问题:你如何衡量你的收益? 𝑥ҧ=𝑥1𝑛1+𝑥2𝑛2+⋯+𝑥𝑘𝑛𝑘 =1σ𝑘 𝑥𝑛 =σ𝑘 𝑛𝑖𝑥 𝑛1+𝑛2+⋯+𝑛𝑘 𝑛 𝑖=1 𝑖 𝑖 𝑖=1 𝑛 𝑖 • 注:𝑛 = 𝑛1 + ⋯ + 𝑛𝑘 • 令𝑤 = 𝑛𝑖，上式可以写成𝑥ҧ = σ𝑘 𝑤 𝑥 𝑖𝑛 𝑖=1𝑖𝑖 描述统计(加权算术平均) • 假定现在你手里有𝑘支股票 • 每只股票有1股 • 第𝑖只股票价格为𝑝 元/股 • 股票收益率依次为:𝑥 ,𝑥 ,⋯,𝑥 12𝑘 • 问题:你如何衡量你的收益? 𝑥ҧ=𝑥1𝑝1+𝑥2𝑝2+⋯+𝑥𝑘𝑝𝑘 =1σ𝑘 𝑥𝑝 =σ𝑘 𝑝𝑖𝑥 𝑝1+𝑝2+⋯+𝑝𝑘 𝑝 𝑖=1 𝑖 𝑖 𝑖=1 𝑝 𝑖 •注:𝑝=𝑝 +⋯+𝑝 𝑖 1𝑘 • 令𝑤 = 𝑝𝑖，上式可以写成𝑥ҧ = σ𝑘 𝑤 𝑥 𝑖𝑝 𝑖=1𝑖𝑖 描述统计(加权算术平均) • 假定现在你手里有𝑘支股票 • 第𝑖只股票，有𝑛𝑖股 • 第𝑖只股票价格为𝑝𝑖元/股 • 股票收益率依次为:𝑥 ,𝑥 ,⋯,𝑥 12𝑘 • 问题:你如何衡量你的收益?  𝑥ҧ = 𝑥1𝑛1𝑝1+𝑥2𝑛2𝑝2+⋯+𝑥𝑘𝑛𝑘𝑝𝑘 = σ𝑘 𝑛1𝑝1+𝑛2𝑝2+⋯+𝑛𝑘𝑝𝑘 𝑤 𝑥 𝑖=1 𝑖 𝑖 •注:𝑤 = 𝑛𝑖𝑝𝑖 ≥0,且 𝑤 +⋯+𝑤 =1 𝑖 𝑛1𝑝1+⋯+𝑛𝑘𝑝𝑘 1 𝑘 描述统计(中位数) • 平均的缺陷: • 被平均了...... • 还没跑赢大盘...... • 对不起，又拖了祖国的后腿...... • 为什么会有这种现象? • 问:平均数的缺陷是什么? • 答:容易受极值影响。 • 问:如何解决? • 答:中位数𝑥𝑀𝐸𝐷(大小排在中间的那个数) • 即𝑥((𝑛+1)/2)或(𝑥(𝑛)+𝑥(𝑛+2))/2 22 补充:分位数 • 分位数: • 假定有数据 𝑥(1),𝑥(2),⋯,𝑥(𝑛) • 中位数:前后各有50%的数 •𝑝分位数:前面有𝑝×100%的数，后面有1−𝑝 ×100%的数 • 上𝑝分位数，下𝑝分位数 • 四分位数(quartile): • Q1:下0.25分位数 • Q3:上0.25分位数 补充:几何均值 • 假定你有一只股票，持有了𝑇天，第𝑖天的收益率为𝑟 ，日均收益率怎么 算? 𝑖 𝑇 • 几何均值(geometric mean/average) • 𝑥1𝑥2⋯𝑥𝑛 1/𝑛 (𝑥𝑖必须为正数) • 适用场合:平均增长率等 • 金融收益率 • GDP 增长率 • 年均增长率 •1+𝑟 1+𝑟 ⋯1+𝑟 =1+𝑟ҧ 12𝑇 补充:基尼系数 • 用于衡量财富/收入分布: 𝐴 𝐴+𝐵 • 取值范围:0~1 • 越靠近0，说明收入差异越小 • 越靠近1，说明收入差异越大 源自wiki 补充:基尼系数 • 右偏数据 • 数据来源:中国家庭金融调查，中国综合社会调查 源自年收入 100 万在目前中国属于什么水平? 补充:恩格尔系数 • 食品支出除以总支出金额 • 根据联合国粮农组织提出的标准，恩格尔系数 在59%以上为贫困，50-59%为温饱，40-50%为小 康，30-40%为富裕，低于30%为最富裕。 • 国家统计局发布《2013年国民经济和社会发展 统计公报》显示，2013年中国农村居民恩格尔 系数为37.7%，城镇居民恩格尔系数为35.0%， 均较上年下降。(大智慧阿思达克通讯社2014 年2月24日讯) 源自wiki 描述统计(中心) ##### Summary Statistics of Audit Data ##### install.packages("prettyR") library("prettyR") Mode(rate.A) order(Audit) sort(Audit) sum(Audit) mean(Audit) median(Audit) min(Audit) max(Audit) quantile(Audit) quantile(Audit,0.25) quantile(Audit,0.75) # randomly assign weights w = runif(length(Audit)) w = w/sum(w) weighted.mean(Audit,w) 单变量数据的描述统计 如何刻画离散(离心) 数据中的信息(变异性) • 极差(range) • 𝑥(𝑛) − 𝑥(1) • 四分位差(inter-quartile range) • Q3 − Q1 源自Statistics for Business and Economics (11th Edition) 数据中的信息(变异性) • 不是K-线图! 源自Statistics for Business and Economics (11th Edition) 数据中的信息(变异性) • 盒子图(箱线图的局限性) 源自网络 数据中的信息(变异性) • 所有的数离中心的距离 • 𝑥1 −𝑥ҧ,𝑥2 −𝑥ҧ,⋯,𝑥𝑛 −𝑥ҧ • 𝑥1 − 𝑥𝑀𝐸𝐷,𝑥2 − 𝑥𝑀𝐸𝐷,⋯,𝑥𝑛 − 𝑥𝑀𝐸𝐷 • 总结所有数据的离心程度 • 1 σ𝑛 |𝑥 − 𝑥∗| (mean of Absolute Deviation) 𝑛 𝑖=1 𝑖 • 1 σ𝑛 𝑥𝑖 − 𝑥∗ 2 (mean of Squared Deviation) 𝑛 𝑖=1 • 𝑥∗可以是𝑥ҧ或𝑥𝑀𝐸𝐷 • 数理性质 • 对任意𝑥，有 1σ𝑛 𝑛 𝑖=1 𝑥𝑖 −𝑥𝑀𝐸𝐷 ≤1σ𝑛 𝑥𝑖 −𝑥 𝑛 𝑖=1 • 对任意𝑥，有 1 σ𝑛 𝑛 𝑖=1 𝑖 𝑛 𝑖=1 𝑖 𝑥 − 𝑥ҧ 2 ≤ 1 σ𝑛 𝑥 − 𝑥 2 数据中的信息(变异性) • 方差(variance) • 𝜎2 = 1 σ𝑛 𝑥 − 𝑥ҧ 2 𝑛 𝑖=1 𝑖 • 标准差(Standard Deviation):𝜎 • 直观解释 • 𝜎 ≥ 0，(𝜎 = 0表示变异程度为0) • 𝜎大小与变异程度大小对应 数据中的信息(变异性) • 过年包饺子的方差 数据中的信息(变异性) • 方差视觉化(上证综合指数2010.1.1-2019.3.20) 数据中的信息(变异性) • 方差的应用 • 在风险管理中的应用(综合考虑收益与波动) • 在质量控制中的应用(生产出来的零件直径) • 在体育训练中的应用 单变量数据的描述统计 分布 正态分布的形状 偏度(skewness) • 左(负)偏与右(正)偏 • 如何刻画偏度: • 𝑆𝑘𝑒𝑤𝑛𝑒𝑠𝑠 = 𝐸 • 具体案例: • 对称分布的偏度均为0 • 常见非对称分布:指数分布等 3 𝑠 𝑥𝑖−𝑥ҧ 峰度(kurtosis) • 高峰分布与低峰分布 •Kurtosis=𝐸 𝑥𝑖−𝑥ҧ 4−3 𝑠 • 注: • 正态分布的峰度为0 • 常见高峰分布:金融数据 高考成绩的分布 https://www.zhihu.com/question/31723932 高考成绩的分布 2017年部分省市高考成绩直方图 数据来源https://www.zhihu.com/question/31723932 高考成绩的分布 数据来源https://www.zhihu.com/question/31723932 正态概率图(QQ图) • 𝑍 ∼ 𝑁(0,1) •𝑋=𝜇 +𝜎𝑍⇒𝑋∼𝑁(𝜇,𝜎2) 𝑋𝑋𝑋𝑋 •𝛼=𝑃𝑍≤𝑐 =𝑃𝜇 +𝜎𝑍≤𝜇 +𝜎𝑐 =𝑃𝑋≤𝜇 +𝜎𝑐 𝛼𝑋𝑋𝑋𝑋𝛼𝑋𝑋𝛼 描述统计(离散及分布信息) ##### Summary Statistics 2 of Audit Data ##### diff(range(rate.A)) diff(quantile(rate.A,c(0.25,0.75))) boxplot(rate.A) boxplot(rate.A,outline=F) var(rate.A) sd(rate.A) co.var<-function(x){ 100*(sd(x,na.rm=TRUE)/mean(x,na.rm=TRUE)) } co.var(rate.A) install.packages("moments") library("moments") skewness(rate.A) kurtosis(rate.A) qqnorm(rate.A,main="Normal Q-Q Plot of rate.A") 数据中的信息(总结) 离散变量 连续变量 集中信息 众数 均值 中位数 变异信息 比例 极差 四分位间距 方差 标准差 形状信息 偏度 峰度 多变量数据的简单分析 表与图 数据中的信息 • 多变量数据存储形式 • (𝑥,𝑦,𝑧,⋯) 或 𝑥,𝑦,𝑧,⋯ 𝑇 源自Statistics for Business and Economics (11th Edition) 数据中的信息(表) • 离散型+离散型(双变量) 源自Statistics for Business and Economics (11th Edition) 变量之间的关系 • 总的位置信息(均值等) • 总的变异信息(极差等) • 刻画(𝑥 , 𝑦 )之间的关系(多元与一元的区别) 𝑖𝑖 • 𝑥 与 𝑦 之间独立(没有关系) • 𝑥 与 𝑦 之间不独立(有关系) • 什么关系? • 如何刻画这些关系? • 如何利用这些关系? • 案例: • 银行核发信用卡(地址+房价) • 精准广告(打标签) • 数学成绩与物理成绩之间的关系 变量之间的关系 • 全班𝑛个人，记下他们的数学与物理成绩 • 第𝑖个人的数学成绩为𝑥 ，物理成绩为𝑦 • 问:总成绩的均值与方差分别是多少? • 均值 1𝑛1𝑛1𝑛 𝑛 ෍ ( 𝑥 𝑖 + 𝑦 𝑖 ) = 𝑛 ෍ 𝑥 𝑖 + 𝑛 ෍ 𝑦 𝑖 = 𝑥 ҧ + 𝑦ത 𝑖=1 𝑖 𝑖 𝑖𝑖 变量之间的关系 • 全班𝑛个人，记下他们的数学与物理成绩 • 第𝑖个人的数学成绩为𝑥 ，物理成绩为𝑦 • 问:总成绩的均值与方差分别是多少? • 方差 𝑛 𝑛 𝑖𝑖 1෍𝑥𝑖+𝑦𝑖− 𝑥ҧ+𝑦ത 2=1෍(𝑥𝑖−𝑥ҧ)+ 𝑦𝑖−𝑦ത 2 𝑛 𝑖=1 𝑛 𝑖=1 1𝑛1𝑛2𝑛 = 𝑛 ෍ 𝑥 𝑖 − 𝑥 ҧ 2 + 𝑛 ෍ 𝑦 𝑖 − 𝑦ത 2 + 𝑛 ෍ ( 𝑥 𝑖 − 𝑥 ҧ ) ( 𝑦 𝑖 − 𝑦ത ) 𝑖=1 𝑖=1 𝑖=1 只与𝑥有关 只与𝑦有关 𝑥与𝑦之间的关系 变量之间的关系 • 协方差(covariance) •𝜎 =1σ𝑛 (𝑥−𝑥ҧ)(𝑦−𝑦ത)(中心化) 𝑥𝑦𝑛𝑖=1𝑖 𝑖 • 相关系数(correlation) 𝜎 • 𝜌 = 𝑥𝑦 (𝜎 与𝜎 分别为𝑋与𝑌的标准差) 𝑥𝑦 𝜎𝜎 𝑥 𝑦 𝑥𝑦 •−1≤𝜌 ≤1 (标准化) 𝑥𝑦 • 性质 • 𝑥 与 𝑦 各自加一个常数，相关系数不变 • 𝑥 与 𝑦 各自乘一个非零常数，相关系数不变 变量之间的关系 • 思考以下情况时 𝑋与𝑌的相关性 • 𝑋与𝑌均为常数(𝜋与𝑒) • 𝑋与𝑌中有一个为常数，另一个不是常数 • 𝑋变化一个单位，𝑌相应变化𝑘个单位(𝑘 ≠ 0) 变量之间的关系 • 备注 • 仅刻画线性相关性，不能刻画非线性相关性 • 例:𝑥=𝑐𝑜𝑠𝜃, 𝑦=𝑠𝑖𝑛𝜃 • 常见图形与对应的Pearson相关系数: 源自wiki 变量之间的关系 • 上证指数前后连续两个交易日涨跌幅度关系 • 相关系数0.0238 变量之间的关系 • 上证指数前后连续两个交易日涨跌额度关系 • 相关系数0.0385 变量之间的关系 • 全球金融市场 变量之间的关系 • 多个变量之间的相关关系:英超球员技术分析 源自网络 变量之间的关系 • 如何解释很重要 • 单店支付金额与商家数量的关系: 变量之间的关系 • 如何解释很重要 • 收入与颜值关系(高跟鞋曲线): 变量之间的关系 • 均值、标准差(方差)、协方差、相关系数等指标的局限性: 变量之间的关系 • 均值、标准差(方差)、协方差、相关系数等指标的局限性: 描述统计(离散及分布信息) ##### covariance and correlation ##### cov(rate.A,rate.B) cor(rate.A,rate.B) cov(SHStock[,4:7]) cor(SHStock[,4:7]) 相关关系与因果关系 相关关系与因果关系 • 注意区分因果关系与相关关系 • 因果=>相关，反之不然
• 相关分为线性相关和非线性相关(本课程内只讨论线性相关)
• 三种类型的因果关系导致的相关关系 • 因为A，所以B，于是有相关关系
• A和B互为因果，于是有相关关系
• A和B都受到C的影响，于是有相关关系

相关关系与因果关系
• 注意区分因果关系与相关关系
源自网络

相关关系与因果关系
• 注意区分因果关系与相关关系
源自网络

多变量分析的复杂性
• 下表中是美国职业棒球联盟2002年到2003年两位球员罗尼·贝利亚德 (Ronnie Belliard)和凯西·布雷克(Casey Blake)的击打率数据。 你觉得谁的表现更好?
年份
2002
运动员
贝利亚德
布雷克
击中数
61
4
击打次数
289
20
平均击打率
0.2111
0.2000
2003
贝利亚德
布雷克
合计 贝利亚德 185 736 0.2514 布雷克 147 577 0.2548
124
143
447
557
0.2774
0.2567

多变量分析的复杂性
• 下表中是美国职业棒球联盟2002年到2003年两位球员罗尼·贝利亚德 (Ronnie Belliard)和凯西·布雷克(Casey Blake)的击打率数据。 你觉得谁的表现更好?
年份
运动员
击中数
击打次数
平均击打率
2002
贝利亚德
61
289
0.2111
布雷克
4
20
0.2000
2003
贝利亚德
124
447
0.2774
布雷克
143
557
0.2567
合计
贝利亚德
185
736
0.2514
布雷克
147
577
0.2548

多变量分析的复杂性
• 录用与性别(辛普森悖论) 某公司男性与女性在过去两年内的录用情况列于下表:
申请人数
录用人数
录用比例
总比例
外勤
内勤
外勤
内勤
外勤
内勤
男
16
4
14
1
87.5%
25%
75%
女
4
16
4
6
100%
37.5%
50%

从认识规律走向利用规律 描述统计->推断统计

描述统计->推断统计
• 统计(收集数据、分析数据、从数据中提取信息)
• 描述统计:整理已有数据资料(结束)
• 概率论:用概率的工具来刻画现有数据的特征(结束) • 推断统计:从已知去探求未知(以后的课程)
• 什么时候才可以从已知推断未知
• 未知(未来)与已知(过去、现在)有联系(如:相关或同分布) • 该联系可以量化表示出来
• 应当有一些规律，不随着是否已知或未知而改变

描述统计->推断统计
http://math.bnu.edu.cn/~lizh/others/yanpoem.htm

描述统计->推断统计
原句
解释
原句
解释
随机
金融收益率
非随意
涨跌比例1:1
人的身高
高个子矮个子相 对都是少数，大 多数人是中等个 子
无序
金融时间序列 (股票价格，期 货价格等)
有序
波动聚集性
红包里面的钱数
毕导定律:每个 红包数服从0.01 到剩下的平均每 人钱数的两倍间 均匀分布。
概率破玄机: 概率论来刻画 随机现象
统计解迷离: 画图，推导， 模拟，验证。

抽样概述

抽样概述
◼ 任意给出一组数据，我们都可以做描述统计，进而说出一系列的故事。 ◼ 但是，为了达到某种目的，存在选择性地挑选对自己有利(或对对手
不利)的数据来说事的可能性。
◼ 为了保证数据分析结果的公正性和可靠性，我们不能简单地对数据采 取“拿来主义”。
◼ 我们需要从源头上对数据的采集、抽样进行“控制”。

抽样概述
HIS $200 BILLION EMPIRE IS BUILT ATOP A MOUNTAIN
OF FAKES.
• 淘宝有假货么?
• 淘宝的假货比例是多少? • 如何抽样调查?

抽样概述
• 案例:工商总局发布2014年下半年网络交易商品定向监测结果。
• 思考:
• 如果你是被检查出有问题的商家，你服不服? • 如果不服，你会怎样质疑工商总局?

抽样概述
• 问题1:你买这些商品想 干嘛?我有很多商品， 为什么只挑一部分?
• 问题2:为什么只挑那一 部分?(工商总局是怎 么抽样的?)
• 是不是刻意知假买假?
• 为什么只抽那么多就下结 论?说不定正好碰上假货 了呢?
• 问题3:从抽得的样本怎 么能推断我所有商品中 假货的比例?你推得准 不准?

抽样概述
• 问题1:
• 问题1.1:抽样的目的是什么? • 问题1.2:为什么要抽样?
• 问题2:抽到的样本是否能代表总体?
• 问题2.1:抽样方式是否合理(或是不合理)?
• 问题2.2:抽样数量是否足够多(或是不够多)?
• 问题3:如何从抽得的样本来估计总体? • 问题3.1:如何估计?
• 问题3.2:估计得准不准?有多准?

抽样概述
• 几个基本概念:
• 个体:收集数据的基本单位
• 总体:所有感兴趣的个体的集合
• 样本:在抽样中被取到的总体的一个子集
• 问:工商总局质检报告中个体、总体、样本都是什么

抽样概述
• 其他日常生活中常见的抽样调查 • PISA调查
• 美国总统选举出口民调(选举)
• 对春晚的评价(好评差评，打分)
• 全国家庭金融调查(开征房地产税的影响)? • 某职位的薪水/行价(人力资源如何定价)? • ……
• 以上抽样调查中个体、总体与样本分别是什么?

抽样调查
问题1.1:抽样的目的是什么?

抽样的目的
• 其他日常生活中常见的抽样调查
• 工商总局的假货调查(均值)
• PISA考试(均值，方差，分布)
• 美国总统选举出口民调(均值，或众数)
• 全国家庭金融调查(均值、方差，分布)?
• 某职位的薪水/行价(均值、中位数、分布)? • ……
• 以上调查的目的都是要计算总体或者总体的某些特征。

抽样的一个附带功能
• 通过抽样，可以直接获取一手数据来验证二手数据
• 获取一手数据的必要性
• 空气质量指数(AQI)是一手数据
• “蓝天”，各种颜色的“预警”是加工后的二手数据
• 任何数据被加工后，一般都会有信息的损失(除非是一一对应的变换)

抽样调查
问题1.2:为什么要抽样?

为什么要抽样
• 为什么要抽样调查来估计总体而不是用全体数据?
• 无法收集到所有总体 • 技术上无法达到
• 成本限制
• 没有必要收集到所有的总体
• 注1:因为不是全部数据，所以通过样本得到的估计往往与真实值有偏差 • 注2:如果样本是随机抽取的，那么通过样本得到的估计往往也是随机的 • 注3:与“大数据”中“全样本”的区别

抽样调查
问题2.1:如何保证抽样准确反映总体 (抽样方式是否合理)

抽样误差
 只要是随机抽样，总会有误差
• 随机抽样过程中，每次抽到的样本是不一样的
• 根据样本算得的结果也往往是不一样的
• 根据样本算得的结果跟总体的真实情况也往往是不一样的
• 这种抽样造成的误差是不可避免的(不论你使用什么方法) • 我们只能尽量以更好的抽样方案来减少抽样误差
 注意区分
• 抽样误差与观测误差

抽样误差
 注意区分抽样误差和抽样过程中的观测误差

抽样方式是否合理
• 两种情况:知道总体信息与不知道总体信息
• 简单随机样本
• 每个样本均以等概率被抽中
• 每个样本是否被抽中与其他样本抽中与否无关
• 总体个数无限情况下
• 有放回和无放回抽样结果不会差异太大
• 总体个数有限情况下
• 有放回与无放回抽样会有差异
• 当总体数量足够多，则可以近似视作无限

抽样方式是否合理
http://www.moe.gov.cn/jyb_xwfb/gzdt_gzdt/s5987/201912/t2 0191204_410707.html

抽样方式是否合理
https://www.oecd.org/pisa/publications/pisa-2018-results.htm

抽样方式是否合理

抽样方式是否合理
• 常见抽样方式(如果我们要在国内做基础教育抽样调查为例)
• 分层随机抽样(stratified random sampling) • 个体依据某一个或几个特征被分成若干层(组)
• 从层(组)里随机抽取
• 减少不确定性(因为层内的样本同质性比较高)
• 整群抽样(cluster sampling) • 先把个体分成若干群
• 一次随机抽取若干群
• 系统抽样(systematic sampling) • 随机编号，随机抽一个
• 以给定间隔往下抽
• 方便抽样(凭便利) • 判断抽样(凭经验)
概率抽样
非概率抽样

抽样方式是否合理
• 时效性(中老年人口)

抽样方式是否合理
• 时效性(淘宝运费险规则的变化)

抽样方式是否合理
• 全面性(抽样不能抽偏样，幸存者偏差) • Wald判断美军飞机的安全防护
• 2018年全国高考理科二卷作文题

抽样方式是否合理
• 全面性(抽样不能抽偏样，幸存者偏差) • Wald判断美军飞机的安全防护
• 2018年全国高考理科二卷作文题

抽样方式是否合理
• 全面性(抽样不能抽偏样，幸存者偏差) • 《柳叶刀》:肥胖的人不容易得老年痴呆
• 《新英格兰医学期刊》:我不这么认为

抽样方式是否合理
• 全面性(抽样不能抽偏样)
• 带鱼侯:《不要因为雾霾而葬送了明天》
事实上，我们把前100个人口大城市空气污染程度进行一个排名的话，你就会 看到:全世界25个污染最严重的城市有13个在印度，兰州是污染最严重的50个 城市里唯一的中国城市;北京的污染程度排名仅为79。——(摘自纽约时报 Gardiner Harris)
今年这位《纽约时报》驻新德里的记者Gardiner Harris放弃驻外记者优厚的 待遇，离开印度回美国工作。临行前他写了一篇文章，详细讲述了自己作这个 决定的原因——新德里的空气污染极度严重，他的小儿子患了严重的哮喘，肺 部功能出现了终生无法恢复的损伤。 他描绘的新德里，让人望而生畏，宛如地 狱。那里一半的学龄儿童肺部受损;孩子们踢球的场边，人工呼吸器扔满了一 地;外籍人士纷纷逃离，甚至让那里的美国学校面临招生不足。他不无悲愤地 这样写道:“在北京，PM2.5值超过500就会登上国际媒体的头条;而德里的数 值是北京的两倍，却基本上无人在意。”

抽样方式是否合理
• 全面性(抽样不能抽偏样)
• 带鱼侯:《不要因为雾霾而葬送了明天》

抽样方式是否合理
• 全面性(抽样不能抽偏样)
• 带鱼侯:《不要因为雾霾而葬送了明天》

抽样方式是否合理
• 全面性(抽样不能抽偏样)
• 带鱼侯:《不要因为雾霾而葬送了明天》

抽样方式是否合理
• 准确性(尤其是主观问题)
• 心理测试中量表设计，多个问题针对一个特征，或是一个问题多次问

抽样方式是否合理
• 有限适用性(科学性) • 普通商品
• 特殊商品

抽样方式是否合理
调查是否反映被调查者的真实意愿(信度)

抽样方式是否合理
• 调查是否反映被调查者的真实意愿(效度)

回顾:区间估计

回顾:区间估计
• 什么是区间估计?为什么要做区间估计?
• 抽样总是有随机性的存在，也就是说，通过抽样得到的估计和真实情况总是有
些不一样
• 很多时候我们不需要知道某个具体的值，我们想知道的是真实的参数(比如正 品率)在哪个范围内
• 区间估计要做的就是估计一下这个“范围”

回顾:区间估计
• 工商总局此次监测共完成92个批次的样品采样，其中有54个批次的样品为 正品，38个批次为非正品。
• 假设工商总局的抽样没有任何问题
• 正品率的估计为0.587
• 有多少把握真实的正品率正好是0.587? • 0!为什么?
• 但是真实的正品率多半(以比较高的概率)在0.587的附近 • 即真值𝜇在区间(0.587 − 𝑑, 0.587 + 𝑑)内的概率比较高
• 即𝑃 0.587−𝑑<𝜇<0.587+𝑑 比较高 • 𝑑越大，区间包含𝜇的概率越大 • 𝑑越小，区间包含𝜇的概率越小 回顾:区间估计 ത • 我们希望估计得到的样本均值𝑋和真实的均值𝜇很接近。 ത • 所谓很接近，从数学上来说，就是给定一个小的正数𝑑，误差 𝑋 − 𝜇 < 𝑑 (𝑑为distance的缩写) തത • 由于𝑋的随机性， 𝑋 − 𝜇 ≥ 𝑑这事总有一定的概率会发生 തത • 但是大多数情况下样本均值𝑋和总体均值𝜇很接近。(即 𝑋 − 𝜇 < 𝑑) ത • 所谓“大多数情况很近”就是指发生的概率比较高，也就是说:𝑃( 𝑋 − 𝜇 < 𝑑)比 较大 • 问:多大算大? • 答:可以事先给定一个概率𝑝，不小于𝑝就算大(𝑝可以取0.9，0.95，0.99，等等) തത • 总结上面的讨论，对于𝑋，有:𝑃 𝑋−𝜇 <𝑑 ≥𝑝 回顾:区间估计 • 以下三种情况等价: ത •𝑋−𝜇<𝑑 ത •𝜇−𝑑<𝑋<𝜇+𝑑 തത •𝑋−𝑑<𝜇<𝑋+𝑑 • 以下三个概率相同: ത •𝑃(𝑋−𝜇 <𝑑) ത • 𝑃(𝜇 − 𝑑 < 𝑋 < 𝜇 + 𝑑) തത • 𝑃(𝑋 − 𝑑 < 𝜇 < 𝑋 + 𝑑) 回顾:区间估计 • 以下说法等价: ഥ • 𝑿与𝝁的距离小于𝒅的概率大于等于𝑝 ത 𝑃𝑋−𝜇<𝑑≥𝑝 • 区间 𝑿−𝒅,𝑿+𝒅 包含𝝁的概率要大于等于𝑝，即 തത 𝑃𝑋−𝑑<𝜇<𝑋+𝑑 ≥𝑝 ഥ • 𝑿落在区间 𝝁−𝒅,𝝁+𝒅 内概率大于等于𝑝，即 ത 𝑃𝜇−𝑑<𝑋<𝜇+𝑑 ≥𝑝 തത • 称𝑝为置信度，(𝑋 − 𝑑, 𝑋 + 𝑑)为置信区间， 𝑑为置信区间半径 • 思考:𝒅与𝒑的关系如何? • 简单起见，我们后面只考虑上述式子中等于的情况 ഥഥ 回顾:区间估计 • 区间估计在解决两个问题: •若𝑝已知，这个区间等于多少?𝑃 𝑋−𝜇 <𝑑 =𝑝，𝑑=? •若𝑑已知，这个概率等于多少?𝑃 𝑋−𝜇 <𝑑 =𝑝，𝑝=? 2 • 由CLT，如果总体均值为𝜇，方差为𝜎 ，则有𝑋 ∼ 𝑁(𝜇, 𝑛 )，即𝜎/ 𝑛 ∼ 𝑁(0,1) തത ത ത ത 2 ത 𝜎 𝑋−𝜇 • 所以𝑃 𝑋 − 𝜇 < 𝑑 = 𝑃 −𝑑 < 𝑋 − 𝜇 < 𝑑 ത =𝑃−𝑑<𝑋−𝜇<𝑑=𝑝 𝜎𝜎𝜎 𝑛𝑛𝑛 回顾:区间估计 标准正态分布 ഥ 𝑃−𝑑<𝑋−𝜇<𝑑=𝑝 𝜎𝜎𝜎 𝑛𝑛𝑛 回顾:区间估计 标准正态分布 ത 𝑃−𝑑<𝑋−𝜇<𝑑=𝑝 𝜎𝜎𝜎 𝑛𝑛𝑛 • 回顾标准正态分布，假如正中间的面积如果是 • 90% • 95% • 99% 对应的 𝑑 分别是多少? 𝜎/ 𝑛 答案: • 1.645(对应90%) • 1.960(对应95%) • 2.576(对应99%) 回顾:区间估计 标准正态分布 ത • 𝑃 − 𝑑 < 𝑋−𝜇 < 𝑑 = 𝑝 • 假定𝑝=0.90，则𝑑 =1.645⇔𝑑=1.645 𝜎 𝜎𝜎𝜎 𝑛𝑛𝑛 𝜎𝑛 𝑛 • 假定𝑝=0.95，则𝑑 =1.960⇔𝑑=1.960 𝜎 𝜎𝑛 𝑛 • 假定𝑝=0.99，则𝑑 =2.576⇔𝑑=2.576 𝜎 𝜎𝑛 𝑛 • 如果𝑛与𝜎已知，则 𝑑 与 𝑝 可以通过上式最后一个等号互推。(一一对应关系) 回顾:区间估计 • 区间半径都有相同的形式𝑑 = 𝑐 • 其中𝑐根据𝑝的不同而取不同的值。 • 以后我们记 • 𝑐𝑁,0.90 = 1.645 • 𝑐𝑁,0.95 = 1.960 • 𝑐𝑁,0.99 = 2.576 • 下标中的𝑁表示正态分布，数字表示置信度 • 有时候简便起见，我们就直接写𝑐，不写下标 𝜎 𝑛 抽样调查 问题2.2:抽样的数量是否足够多(或不够多) 抽多少? 抽多少?  要回答“抽多少”需要考虑的因素:  抽样的目的是什么?  估计  手段为目的服务  你对抽样后估计的精度要求有多高?  要求越高，样本量越大  你的预算有多少?  巧妇难为无米之炊  待抽样的总体长什么样?  如果大家都长一样，抽一个就够了 抽多少?  问题与反问题:  给定估计精度至少需要多少样本?抽这么多需要多少成本?  给定成本条件下，最多能抽多少?能达到什么样的精度?  之所以“抽多少”这个问题不好回答，是因为:  不知道总体的分布情况  不知道对估计的精度要求是什么  不知道预算是多少?  此外，由于抽样有随机性，所以抽出来的结果也有随机性，必须考虑到这种随机性的影 响。(比如不同的抽样方式会有不同的随机结果)  既然最终目的是估计，我们先放下抽样，看看估计。 抽多少? • 回顾一下区间估计 大样本(𝒏 > 𝟑𝟎)
小样本(𝒏 ≤ 𝟑𝟎)
𝜎已知
(𝑥ҧ − 𝑐𝑁,𝑝 × 𝜎 , 𝑥ҧ + 𝑐𝑁,𝑝 × 𝜎 ) 𝑛𝑛
总体分布是正态时
(𝑥ҧ − 𝑐𝑁,𝑝 × 𝜎 , 𝑥ҧ + 𝑐𝑁,𝑝 × 𝜎 ) 𝑛𝑛
总体分布不是正态时 很复杂
𝜎未知
(𝑥ҧ−𝑐𝑁,𝑝× 𝑠 ,𝑥ҧ+𝑐𝑁,𝑝× 𝑠) 𝑛𝑛
总体分布是正态时
(𝑥ҧ − 𝑐𝑡(𝑛−1),𝑝 × 𝑠 ,𝑥ҧ + 𝑐𝑡(𝑛−1),𝑝 × 𝑠 ) 𝑛𝑛
总体分布不是正态时 很复杂

抽多少?
• 考虑以下三个调查案例:
• 调查1:调查某个职位的人员平均薪资(连续)
• 调查2:调查支持开征房产税的人比例(离散)
• 调查3:调查准备生二胎的适龄夫妻比例(离散)
• 调查的范围在以下两个城市: • 上海市，常住人口为2000万
• 某中等城市，常住人口为200万
• 请问:
• 你认为应该调查多少人?
• 如果其他条件都一样，上海市要调查10000人，则那个中等城市该调查多少人?

抽多少?
• 因为抽样是为了估计总体的某个参数(比如均值)，所以抽样多少与总 体有关。
• 由于抽样有随机性，所以算出来的结果也有随机性，必须考虑到这种随 机性的影响。
• 抽多少?
• 与要求的精度有关系
• 与总体的离散程度有关系
• 讨论:与总体的大小有没有关系?

抽多少?
തത
⚫区间 𝑋−𝑑,𝑋+𝑑 包含𝜇的概率要大于等于𝑝，即
തത ⚫𝑃𝑋−𝑑<𝜇<𝑋+𝑑 ≥𝑝 ⚫ 我们可以把置信区间半径𝑑看做是精度 ⚫ 包含𝜇的概率𝑝为置信度 ⚫ 实际工作中我们既要求精度高，也要求置信度高。 ⚫ 比如，如果我们希望满足下面两个条件，需要的样本量至少要多少? ⚫ 精度控制在$1000之内(𝑑 ≤ 1000) ⚫ 置信度要达到99%(𝑝 > 0.99)

抽多少?
⚫ 我们希望样本量越多越好
⚫ 增加样本量的好处:
⚫ 直观解释:样本越多，结果越精确，越可信。
ത𝜎ത𝜎 ⚫ 数理解释: 𝑋−𝑐 𝑛,𝑋+𝑐 𝑛
⚫ 保持置信度(𝑐)不变，𝑛越大，上述区间半径就越小(越精确) ⚫ 保持精度(半径𝑐 𝜎 )不变， 𝑛越大，则𝑐也越大，置信度就越高
𝑛
⚫ 但是样本不能无限制地多，因为收集样本是要成本的，我们总希望节约成本。
⚫ 我们希望抽取的样本数量在满足对于精度和置信度的要求的基础上，越少越好。

抽多少?
⚫ 给定精度𝑑，置信度𝑝，则至少需要多少样本?
⚫ 大样本，总体标准差已知时
2 𝑁,𝑝 𝑛 𝑑
⚫ 𝑑≥𝑐 𝜎 ⇔𝑛≥ 𝜎𝑐𝑁,𝑝
⚫ 大样本，总体标准差未知时
2 𝑁,𝑝 𝑛 𝑑
⚫ 𝑑≥𝑐 𝑠 ⇔𝑛≥ 𝑠𝑐𝑁,𝑝

抽多少?
⚫ 给定精度𝑑，置信度𝑝，则至少需要多少样本?
⚫ 小样本，总体是正态分布，标准差已知时
2 𝑁,𝑝 𝑛 𝑑
⚫ 𝑑≥𝑐 𝜎 ⇔𝑛≥ 𝜎𝑐𝑁,𝑝
⚫ 小样本，总体是正态分布，标准差未知时
2 𝑡(𝑛−1),𝑝 𝑛 𝑑
⚫ 𝑑≥𝑐 𝑠 ⇔𝑛≥ 𝑠𝑐𝑡(𝑛−1),𝑝

抽多少?
◼例子: EAI公司抽样问题
◼ EAI公司准备从市场上招一名管理人员，想知道目前就业市场上该岗位的平均工
资是多少?
◼ 人事经理决定随机抽取30名类似的管理人员，组成一个随机样本
◼ 抽样进行前，并不确定哪30名管理人员将被抽到，故他们的年薪值被视为随机 变量，记为𝑋 ,⋯,𝑋
1 30
◼ 抽样进行后，有具体的30名管理人员将抽到，他们的年薪值被记录下来，成为 30个具体数值 𝑥1,⋯,𝑥30，并可以记录在数据表中

抽多少?
⚫例子: EAI公司抽样问题 ⚫ 30位管理人员的工资
⚫ 样本均值51814，样本标准差3348

抽多少?
⚫ 精度控制在$1000以内，置信度要达到99%，需要多少的样本量? ⚫ 计算:
⚫ 𝑐 𝑠 =𝑑<1000⇔𝑛≥ 𝑐𝑠 2 = 2.58×3348 2 ≈74.61 𝑛 1000 1000 ⚫ 显然，我们至少要75个样本 抽多少? ⚫ 屠呦呦工作的有关新闻报道: ⚫ 经研究发现，双氢青蒿素在免疫领域具良好的双向调节作用，既能降低B细胞高反应性以减少免疫复 合物沉积所致的自身免疫病，又可提高T细胞的免疫功能。在北医有关部门支持下，我们已将双氢青 蒿素用于治疗红斑狼疮和光敏性疾病。现已获国家食品药品监督管理局的“药物临床研究批件” (2004L02089)和中国发明专利(专利号:ZL 99103346.9)。经临床100例疗效初步观察，总有效 率94%，显效率44%。 抽多少? ⚫ 例:双氢青蒿素的免疫试验，临床100例疗效初步观察，总有效率94%， 显效率44%。 ⚫问: ⚫ 双氢青蒿素的显效率的点估计值为多少? ⚫ 为双氢青蒿素的显效率构造一个置信水平为95%的置信区间。 抽多少? ⚫ 分析思路: ⚫ 假定显效率的比率为𝑞，随机选取一个样本。如果第𝑖个样本呈显效率， 则记为𝑋 = 1;否则，记为𝑋 = 0。显然，𝑋 服从参数为𝑞的Bernoulli 𝑖𝑖𝑖 分布。 ⚫ 视病人为独立同分布随机样本 ⚫ 显效率的点估计:𝑞ො = 𝑥ҧ = 44 (大数定律) 100 ⚫ 注:此处用𝑞是因为𝑝已经被用来代表置信度了 抽多少? ⚫ 假定显效率为𝑞，随机抽取1个样本。如果第𝑖个病例呈显效率，则记为𝑋𝑖 = 1;否则， 记为 𝑋𝑖 = 0。显然，𝑋𝑖服从参数为𝑞的Bernoulli分布。 ⚫ 视病人为独立同分布随机样本 ⚫随机做了𝑛个样本，当𝑛足够大(不小于30)时，根据C.L.T.，有𝑋∼𝑁 𝑞, ⚫ 置信水平为95%的置信区间: ⚫ ( 𝑞ො − 1 . 9 6 𝑞ො 1 − 𝑞ො , 𝑞ො + 1 . 9 6 𝑞ො ( 1 − 𝑞ො ) ) (0.3427, 0.5373) = (34.27%, 53.73%) 如果嫌宽? 如何改进? ത 𝑞 1−𝑞 𝑛 𝑛𝑛 抽多少? 𝑞(1−𝑞) 𝑛 ⚫ 置信区间半径:𝑑 = 𝑐 ⚫ 如果我们要求置信区间半径𝑑要小于某个上限𝐷 ⚫𝑐 𝑞(1−𝑞)=𝑑≤𝐷⇔𝑛≥𝑐2𝑞 1−𝑞 𝑛 𝐷2 ⚫问:真实的 𝑞 不知道，如何处理? ⚫答: ⚫ 方案1:预抽样去估算𝑞 ⚫ 方案2:不等号右边最多是多少?只要比最多都多，那就是足够多。所以令 𝑞 = 0.5，这样可使得𝑛足够多，可以满足任意0 ≤ 𝑞 ≤ 1的情况。 抽多少? ⚫ 假定屠呦呦希望构造一个95%的置信区间，控制实际误差不超过0.025，应该做多少 病例(样本)? ⚫ 解: ⚫ 𝑐 = 1.96 (置信度95%) ⚫ 𝐷 = 0.025 ⚫𝑐 𝑞(1−𝑞)≤𝐷⇔𝑛≥𝑐2𝑞1−𝑞 𝑛 𝐷2 ⚫𝑛≥1.962×0.5× 1−0.5 =1536.64≥𝑐2𝑞 1−𝑞 0.0252 𝐷2 (0≤𝑞≤1) ⚫ 思考:如果要求误差不超过0.0025呢?(精度增加一个数量级) ⚫ 答:样本量为原来的100倍(样本量增加两个数量级) 抽多少? • 影响样本量的因素: •𝑛≥ 𝑐𝜎 2 𝐷 • 𝑛 ≥ 𝑐2𝑝 1−𝑝 𝐷2 (连续情况) (比率情况) • c:一个与置信度对应的常数 • 𝐷:要求的精度 • 𝜎:总体的(或样本)标准差 • 注意:没有涉及到总体的大小 𝑁! • 假定有一个2000万人的总体和200万人的总体 • 要求的精度、置信度一样，总体的标准差也一样 • 那么需要调查的样本量也一样!(跟总体大小没有关系!) 补充:假设检验 假设检验是什么? • 在英国剑桥一个夏日的午后，一群大学的绅士和他们的夫人们，还有来访者， 正围坐在户外的桌旁，享用着下午茶。在品茶过程中，一位女士坚称:把茶加 进奶里，或把奶加进茶里，不同的做法，会使茶的味道品起来不同。在场的一 帮科学精英们，对这位女士的“胡言乱语”嗤之以鼻。这怎么可能呢?他们不 能想象，仅仅因为加茶加奶的先后顺序不同，茶就会发生不同的化学反应。然 而，在座的一个身材矮小、戴着厚眼镜、下巴上蓄着的短尖髯开始变灰的先生， 却不这么看，他对这个问题很感兴趣。他兴奋地说道:“让我们来检验这个命 题吧!”并开始策划一个实验。在实验中，坚持茶有不同味道的那位女士被奉 上一连串的已经调制好的茶，其中，有的是先加茶后加奶制成的，有的则是先 加奶后加茶制成的。 • 接下来，在场的许多人都热心地加入到实验中来。几分钟内，他们在那位女士 看不见的地方调制出不同类型的茶来。最后，在决战来临的气氛中，蓄短胡须 的先生为那位先生为那位女士奉上第一杯茶，女士品了一小会儿，然后断言这 一杯是先倒的茶后加的奶。 这位先生不加评论地记下了女士的说法，然后，又 奉上了第二杯...... https://baike.baidu.com/item/%E5%A5%B3%E5%A3%AB%E5% 93%81%E8%8C%B6/2921681?fr=aladdin 假设检验是什么? 假定一共试验了10杯 • 问:假定一位女士声称她具有区分能力，而且答对了全部10杯，那么我们是否 该相信她? • 答:我们倾向于相信她很善于区分，或者说不能排除她善于区分的可能性。 • 问:假定一位女士声称她没有区分能力，但是她答对了全部10杯，那么我们是 否该相信她? 10 • 答:因为如果她没有区分能力，纯靠瞎猜的话，猜中全部10杯的概率是1/2 。 我们觉得她很可能在谦(sa)虚(huang)。她所谓的“没有区分能力”多半是 有问题的。 假设检验是什么? • 假设检验的基本思想:反证法! • 假设有一个论断(命题)，我们管它叫A。(事先我们并不知道A是对的 还是错的。) • 我们想验证A是不是对的。 (用学术的语言来说，真的) • 我们暂且假定A是对的，在其基础上做进一步的演绎和推导，得到结论B。 • 我们发现，实际上B成立的概率非常小。 • 这时候，我们会怀疑:命题A会不会其实是错的? 检验的种类 ⚫ 对单样本参数的检验: ⚫ 检验均值是否等于某个特定的数(不同备择假设) ⚫ 对双样本参数的检验: ⚫ 检验均值是否相同 ⚫ 对分布的检验: ⚫ 检验数据是否服从某个特定的分布(比如正态) ⚫ 检验两组数据是否服从相同分布 检验的种类 K-S检验 ⚫ 如果两个分布的分布函数差的最大值太大，则这两个分布不属于同一 个分布。