理工信（上）开卷代码

写在开头

本笔记以非博雅班的课程内容为主体，理论与代码内容大多来自PPT，一些内容来自ai生成。
代码和理论部分感谢踏实勤劳又肯干的kimi，忙前忙后为我讲解生统部分的基础概念，还为我提供例子，进行代码运行与计算，感谢他的不离不弃；同时代码部分也需要感谢deepseek，虽然我也忘了哪些部分是他写的；也感谢不登录不给问的豆包，谢谢他的陪伴与支持。

概率论基础

事件与运算

PPT1的8-30页

样本空间 (sample)：随机试验所有可能的结果构成的集合
事件(event)：样本空间 S 的任一子集 E
- 基本事件：在一定范围内不能再分解的事件
- 复杂事件：由基本事件组成的事件
- 必然事件：一定条件下必然发生的事件，记为 Ω
- 不可能事件：一定条件下肯定不发生的事件，记为 Ø
事件关系
韦恩图 (Venn diagram)是一种用来阐述事件之间的逻辑关系的非常有效的几何表示方法
- 包含关系：若事件A发生必然导致事件B发生，则称事件A包含于B，或事件B包含 A，记作B $\supset$ A或A $\subset$ B
- 等价关系：若事件A包含于事件B，同时，事件B包含于事件 A，即B $\supset$ A或A $\subset$ B，则称事件A与B等价，记作A=B
  A=B只是说A以及B两个随机事件同时发生或者同时不发生，并不是说两个事件是同一事件
- 事件的并 (和)：若事件C表示A、B中至少有一个发生
  则称事件C为A与B的并事件
- 事件的差：若事件表示B发生但事件A不发生，则称该事件为B与A的差，记为B−A
- 事件的交：若事件C表示A与B同时发生，则称事件C为 A与B的交事件
- 互不相容事件：若事件A和事件B不能同时发生，则称事件A与B为互不相容事件，又称为互斥事件
  如果 n 个互斥事件的并为必然事件，那么称这n个事件构成互斥完备群
- 对立事件：若某事件表示事件A必然不能发生，则称此事件为A事件的对立事件(或事件A的补)
  对立关系：若事件A与B有一个发生且仅有一个发生，则称事件A与B互为对立事件
- 独立事件：事件B的发生不影响事件A的发生的可能性P(AB) = P(A)P(B)
  - 应用
    - 数据是否独立：影响统计方法的选择
      - 成组设计：试验对象随机地分配到处理组和对照组，两组数据是独立的；
      - 配对设计：试验对象在处理前、后分别进行了测量，测量前和测量后的数据不独立。
    - 变量间是否相关(独立则不相关)
      - 发病率与性别是否相关；痊愈率与治疗方法是否相关
事件运算的性质

PPT1的26-30页

德摩根定律DeMorgan’s laws：$\overline{AB}$ = $\overline{A}$ + $\overline{B}$

概率

PPT1的31-41

概率(probability)：刻画事件发生可能性大小的数量指标
古典概率模型：总数为n的等概率事件完备组中，若事件A包含其中的m个等概率基本事件，则事件A的概率为$m/n$
- 条件：完备性、互不相容性、等可能性
- 缺陷：首先需要验证“互不相容”和“等可能性”的条件，需要看基本事件的个数是否是有限的
频率 (frequency)：通过大量的重复试验来探寻不同随机事件发生的可能性大小
设 n 次试验中事件 A 发生了 m 次，则比值 m /n 称为
随机事件A发生的频率
在一个固定的数值附近摆动，呈现出一定的稳定性，是概率的近似
两者关系

PPT1的38页
概率的公理化定义：非负性、正规性、可加性

PPT1的39页
概率的基本性质

PPT1的40页
概率的运算
- 加法公式
  若事件A和B为互不相容的随机事件,则有P(A+B)=P(A)+P(B)
  对于任意的两个随机事件A和B，有P(A+B)=P(A)+P(B)-P(AB)
- 乘法公式：P(AB) = P(B)P(A|B)
- 条件概率：度量两事件的关联性
  在已知 A 事件发生的条件下, B 事件发生的概率记作 P(B|A)，称为在A发生的条件下，B发生的条件概率
  P(B|A)=P(AB)/P(A) P(A)≠0
  - AB互不相容：P(A|B) = 0
  - B包含于A：P(A|B) = 1
  - 医学中的条件概率
    - A = 检查阳性，B = 确诊
      阳性预测值 P+ = P(B|A)
      阴性预测值 P- = P($\bar{B}|\bar{A}$)
    - 灵敏度：发生疾病伴随的症状或检验阳性的概率P(A|B)
    - 特异度：无病时不出现症状或检验阴性的概率P($\bar{A}|\bar{B}$)
- 全概率公式：P(B) = P(A1)P(B|A1)+P(A2)P(B|A2)
  - 全概率 P(B) 被分解成了许多部分之和
  - 从另一个角度看把Ai看做导致事件B发生的一种途径。对不同途径 Ai，发生的概率即条件概率P(B|Ai)各不相同，而采取哪种途径却是随机的。P(B)等于诸P(B|Ai)以P(Ai)为权的加权平均
- 贝叶斯公式：$P(A|B) = \frac{P(A)P(B|A)}{P(B)}$
  或者写作$P(Aj|B) = \frac{P(Aj)P(B|Aj)}{\sum{P(Ai)P(B|Ai)}}$
  $后验概率P(Ai|B)=\frac{先验P(Ai)*似然P(B|Ai)}{边际概率P(B)}$
  +案例：肺癌诊断决策
  
  PPT2.1的32-36页

随机变量

概念及其分布

定义：设E是随机试验,其样本空间是Ω={ω}，如果对于每一个ω∈Ω有一个实数X(ω)和它对应，那么得到一个定义在Ω上的单值实函数X(ω)，称X为定义在上的随机变量
- 当随机事件为定性事件时，随机变量取值仅起到编码识别作用
实质：用确定且唯一的数字表示各试验结果
- 随机变量的定义域不一定是数集
- 随机变量取特定值的概率与对应的随机事件的概率一致
与普通函数的区别
- 随机变量随着试验的结果不同而取不同的值，因此在试验之前只能知道它可能的取值范围，不能预先知道它取什么值，但是试验的各结果的出现有一定的概率，所以随机变量取值也有一定的概率
- 普通函数是定义在实数轴上的，而随机变量是定义在样本空间上的
离散型随机变量
- 若随机变量X存在有限或无限可列个可能的取值，则称X为离散型随机变量，P(X=xk) = pk，且总和为1
- 用表格表示随机变量→概率分布表/分布列

PPT3的9-10页

分布

PPT3的11-52页
- 两点分布（01分布）
X 0 1

P 1-p p
- 伯努利实验：在相同试验条件下相互独立且重复进行的一种随机试验，该随机试验只有两种可能结果，发生or不发生。假设该项试验独立重复地进行了n次，那么就称这一系列重复独立的随机试验为n重伯努利试验
- 二项分布：$X \sim B(n,p)$
  在 n 次独立试验中事件 A 发生 k 次的概率：$P(X=k) = C_{n}^{k}p^k(1-p)^{n-k}$
  k = 0,1,……,n 0<p<1
  - 与两点分布关系
    设随机变量X1,X2,……,Xn服从参数为p的两点分布且相互独立，则$X = \sum_{k=1}^{n}{X_k} \sim B(n,p)$
- 泊松分布：$X \sim P(λ)$或$X \sim Pois(λ)$或$X \sim \pi(λ)$
PPT3的24-27页

$P(X=k) = \frac{\lambda^ke^{-\lambda}}{k!}$
k = 0,1,2,…… λ>0
泊松分布不仅可在试验数N很大时作为二项分布的近似分布，还可以用于描述稀疏事件的概率
例如: 医院服务台在单位时间内的总呼叫数
- poisson定理：在n重伯努利试验中，事件A在一次试验中出现的概率为$p_n$(与试验总数n有关),如果当n趋于无穷时，$np_n$趋于常数λ，则有：
  $lim_{n\to \infty}C_{n}^{k}p_n^k(1-p_n)^{n-k} = \frac{\lambda^ke^{-\lambda}}{k!}$
- 几何分布（不考）：$X \sim G(p)$
PPT3的18页
- 超几何分布（不考）：$X \sim H(N,M,n)$
PPT3的19-21页
- 负二项分布（不考）
分布函数
假设 X 是一个随机变量，t是任意的实数，则函数F(t)=P(X≤t)称为随机变量X的分布函数
随机变量的分布函数与概率分布是不损失信息的等价变换,分布函数将所有离散取值的分析问题改变为有限个区域中的分析问题
- 性质
  - F(x)是一个单调不减的右连续函数：F(x+)=F(x)
  - 0≤F(x)≤1，且有$F(-\infty) = lim_{x\to-\infty}F(x) = 0$$F(+\infty) = lim_{x\to+\infty}F(x) = 1$
  - $P(x_1<X≤x_2) = F(x_2)-F(x_1)$
概率分布

PPT3的35-45页

X	0	1
P	1-p	p

设随机变量 X 的分布函数为 F(t)，若存在非负可积函数 f (x)，使得$F(t) = P(X≤t) = \int_{-\infty}^{t}f(s)ds$，则称X为连续性随机变量，f(x)称为随机变量X的概率密度函数

性质
- $f(x)≥0 \forall{x\in (-\infty,+\infty)}$
- $\int_{-\infty}^{+\infty}{f(x)dx} = 1$
- $P(x_1<X≤x_2) = \int_{x_1}^{x_2}f(x)dx$
- $F’(x) = f(x)$
分类
- 均匀分布：$X \sim Unif(a,b)$
  若连续型随机变量 X 的概率密度为$f(x) = \frac{1}{(b-a)} ,a≤x≤b$，则称 X 服从区间 [a, b]上的均匀分布(Uniform distribution)
- 正态分布：X $\sim$ N($\mu,\sigma^2$)
  若连续型随机变量 X 的密度为$f(x) = \frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x-\mu)^2}{2\sigma^2}},\sigma>0$，则称 X服从参数为𝝁,$\sigma^2$ 的正态分布(NormalDistribution)
  - 单峰对称：正态分布有一个对称中心，且取值偏离中心值的概率与距离的平方成反比
  - 标准正态分布：N(0,1)
    $\phi(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{x}},-\infty<x<+\infty$
  - 正态分布的标准化：设随机变量$X \sim N(𝝁,𝝈^𝟐)$ ，则经过标准化变换$Y=\frac{X-\mu}{\sigma}$可得到$Y \sim N(0,1)$
  - 3σ原则：$P(\mu-3\sigma≤X≤\mu+3\sigma) = 0.9974$
- 其他自学内容
PPT3的46-51页

数字特征

期望E(X)

PPT4.2的6-16页

设离散型随机变量 X 的分布列为$P(X = x_k) = p_k, k = 1,2,……$，若级数$\sum_{k=1}^{\infty}{x_kp_k}$绝对收敛则称该级数为X的数学期望

本质：本质是加权平均，几何意义是重心
连续型的期望：$E(X)=\int_{-\infty}^{+\infty}xf(x)dx$
性质
- 常数期望为0
- E(CX) = CE(X)
- E(X+Y) = E(X) + E(Y)
- E(kX+bY) = kE(X) + bE(Y)
- 相互独立变量：$E(X1*X2) = E(X1)*E(X2)$
  +不同分布的期望

PPT4.2的12-15页

+ 两点分布：EX = p
+ 二项分布：EX = np
+ 泊松分布：EX = λ
+ 正态分布：EX = μ

方差：随机变量相对于中心值的聚散程度

PPT4.2的17-23页

设 X 是一个随机变量， X 的数学期望存在E(X) = μ，若 $E(X‒μ)^2$ 存在且为有限值，则称其为X的方差，记作D(X) 或 Var(X)，$DX = \sigma^2 =E[X-E(X)]^2$

离散型：$DX = \sum{(X_k-E(X))^2}p_k$
连续型：$DX = \int_{-\infty}^{+\infty}(x - E(X))^2f(x)dx$
常用计算公式：$D(X) = E(X^2)-E^2(X)$
性质
- 常数方差为0
- $D(CX) = C^2D(X)$
- 相互独立变量：D(X+Y) = D(X)+D(Y)
标准差：保持量纲一致
$S_d = \sqrt{DX}$
变异系数：比较不同量纲的随机变量的离散程度
变异系数可以消除单位或平均数不同对两个或多个指标变异程度比较的影响
变异系数是衡量数据离散程度的相对指标，常用于比较
不同数据集或不同单位的变异程度
$CV = \frac{\sqrt{DX}}{EX}$
- 红细胞分布宽度RDW：(红细胞平均体积的标准差/平均红细胞体积 MCV ) × 100%
  它是血液计数报告中一项用来测量红细胞变量宽度的参数。如遇到一些变异的情况，当然细胞体积的变量将会明显受到影响。 RDW值越高表明体积的变异就越大
分位数

PPT4.2的27-32页

上侧分位数
对于任意类型的随机变量X，如果能找到数t，使得$P(X≥t) = \alpha$，则称t为变量X的上侧𝜶分位数，记作$x_𝜶$，即 t = $x_𝜶$
双侧分位数
对于任意类型的随机变量X，如果能找到数t，使得$P(|X|≥t) = \alpha$，则称t为随机变量X的双侧𝜶分位数，记作$t = x_{\alpha/2}$
离散随机变量函数的分布

PPT4.2的33-37页
随机变量和的分布

PPT4.2的38-42页

二项分布：$Z \sim B(n1+n2,p)$
泊松分布：$Y \sim Pois(λ1+λ2)$
正态分布：$Z \sim N($\mu_1+\mu_2,\sigma_1^2+\sigma_2^2$)$

随机向量
- 二维随机向量：以随机变量X,Y作为分量所构成的向量(X, Y)
- 联合分布函数：设(X,Y)为二维随机向量，对任意实数 x, y, 二元函数F(x,y)=P(X≤x,Y≤y)，称为二维随机向量的联合分布函数
PPT4.2的44-45页：连续性和离散型的联合密度函数
随机变量的独立性

PPT4.2的46页

设X和Y为两个随机变量，若对任意 a<b,c<d,事件{a<X<b} 与事件{c<Y<d}均相互独立，则称随机变量X和Y相互独立

定理
设 (X, Y)为两维随机向量, X 和 Y 相互独立
- (X, Y)为离散型随机向量
  P(X = xi,Y = Yj) = P(X=xi)P(Y=yj)
- (X, Y)为连续型随机向量
  f(x,y)=fX(x)fY(y)
总体协方差
- $\mu_x = EX,\mu_y = EY$
- $Cov(X,Y) = E[(X-\mu_x)(Y-\mu_y)]$
- $Cov(X,X) = E[(X-\mu_x)(X-\mu_x)] = DX$

大数定理

PPT4.2的48-52页

贝努利大数定理：设m是n次独立重复试验中事件A发生的次数， p是每次实验A发生的概率，则对任意正数ε,总有$lim_{n\to\infty}P(|\frac{m}{n}-p|<\epsilon) = 1$
切比雪夫大数定理：设随机变量Xi相互独立且具有相同的数学期望μ和方差σ^2，若令$\bar{X} = \frac{1}{n}\sum_{i=1}^n{X_i}$，则对任意整数$\epsilon$，总有：$lim_{n\to\infty}P(|\bar{X}-\mu|<\epsilon) = 1$
- 意义：n足够大的时候，算术平均值为一常数，接近与期望

中心极限定理

PPT4.2的53-60页

假设被研究的对象可由大量独立因素的总和解释，其中每个因素对于总和的作用都很微小且可比较，则当因素数目很大时就可以认为这个对象服从正态分布的
如果随机变量序列$X_k(k=1,2,……)$独立同分布，并且有$E(X_k) = \mu$，$D(X_k) = \sigma^2$，则有：$lim_{n\to\infty}P(\frac{\bar{X}-\mu}{\frac{\sigma}{\sqrt{n}}}) = \int_{-\infty}^{t}{\frac{1}{\sqrt{2\pi}}e^{-\frac{u^2}{2}}du}$

一般形式：$lim_{n\to\infty}P(\frac{f(X_n)-E[f(X_n)]}{\sqrt{D[f(X_n)]}}≤t) = \Phi(t)$
二项分布的正态近似
$lim_{n\to\infty}P(\frac{X-np}{\sqrt{npq}}) = \int_{-\infty}^{t}{\frac{1}{\sqrt{2\pi}}e^{-\frac{u^2}{2}}du}$
泊松分布的正态近似
$lim_{n\to\infty}P(\frac{X_n-n\lambda}{\sqrt{n\lambda}}) = \int_{-\infty}^{t}{\frac{1}{\sqrt{2\pi}}e^{-\frac{u^2}{2}}du}$
应用
- 当样本量比较大时，可以用来对给定的随机变量的分布进行近似
- 当研究对象可以表示为独立同分的总和时，可以用正态分布来近似计算

数据

数据的类型与整理

PPT5的5-8页

定义：数据data是对客观事物计量的结果
形式：数字、图像、文本、音频、视频
数据是统计分析的基础
类型
- 定性数据qualitative data
  - 定类数据categorical data：性别、血型
    定义: 事物的属性数据，多表现为文字类型的无序类别数据
    特点: 无大小关系不能排序，无加减运算
  - 定序数据ordinal data：疗效、病情
    定义：是对事物之间等级或顺序差别的观测数据，
    数据表现为有序类别
    特点：有顺序、大小关系，但不能加减运算
- 定量数据quantitative data：年龄、身高、体重
  定义: 对事物按度量单位进行计量或观测的数据，
  表现为具体的数值
  特点: 有数值大小关系，能排序和加减运算
  - 离散型变量：观测只取有限个或可列个值
  - 连续型变量：取无穷多个值，取值连续不断不能一一列举
数据的统计整理和图示

PPT5的9-21页

数据的统计整理就是根据统计研究的任务，对搜集到的数据资料进行科学的汇总和处理，使数据资料系统化，以反映研究总体的特征、规律和趋势

整理步骤：审核并清洗数据→分类/分组→汇总→统计图表报告
数据统计整理
- 定性数据：分类整理
- 定性资料：常用指标
  频率、构成比、率和比
- 定量数据：分组整理
  目的：展示定量数据的分布规律、特征和分布类型
  步骤：确定组数→确定组距→计算频数→频数分布表→统计图示
  - 组数：k = 1 +lnN/ln2
  - 组距：d = [最大值-最小值/组数]（取整）/分位数
定量数据分布特征的统计描述
分布的中心、散布程度、形状：见【数据的概括性度量】部分

数据概括性度量

统计量意义

PPT5的22-28页

中心趋势：均值、中位数
散布程度：方差、标准差、变异系数、极差、四分位间距
形态：偏度、峰度
均值、方差与极值、极差

1 2	`data = c(1,2,3,4) mean(data)`

1	`## [1] 2.5`

1	`exp(mean(log(data)))# 几何平均`

1	`## [1] 2.213364`

1	`var(data)# 方差 1/(n-1)`

1	`## [1] 1.666667`

1	`sd(data)# 标准差 1/(n-1)`

1	`## [1] 1.290994`

1 2	`se = sd(data)/sqrt(length(data))# 是sqrt(D(X_bar)) max(data) - min(data)`

1	`## [1] 3`

1	`range(data)`

1	`## [1] 1 4`

原点矩

1 2	`data = c(1,2,3,4) library(moments)`

1	`## Error in library(moments): 不存在叫'moments'这个名称的程序包`

1	`moment(data)`

1	`## Error in moment(data): 没有"moment"这个函数`

频数分布

data = c(1,2,3,4)
breaks = seq(0, 5, length.out = 3)# 区间数是length.out - 1
data.cut = cut(data, breaks, right=FALSE)
data.freq = table(data.cut)
data.freq

1
2
3

## data.cut
## [0,2.5) [2.5,5) 
##       2       2

分位数
设随机变量X，给定α(0<α<1)，如有数xα，使得$P(X≥x_\alpha)≥\alpha$$P(X≤x_\alpha)≤1-\alpha$，则称xα为随机变量X的上侧 α 分位数

1 2	`data = c(1,2,3,4) quantile(data)`

1 2	`## 0% 25% 50% 75% 100% ## 1.00 1.75 2.50 3.25 4.00`

1	`median(data)`

1	`## [1] 2.5`

1	`quantile(data,0.5)`

1 2	`## 50% ## 2.5`

1	`IQR(data)# 四分位间距`

1	`## [1] 1.5`

变异系数

1 2	`data = c(1,2,3,4) sd(data)/abs(mean(data))`

1	`## [1] 0.5163978`

偏度与峰度

1	`library(moments)`

1	`## Error in library(moments): 不存在叫'moments'这个名称的程序包`

1	`skewness(iris$Sepal.Length)`

1	`## Error in skewness(iris$Sepal.Length): 没有"skewness"这个函数`

1	`kurtosis(iris$Sepal.Length)-3`

1	`## Error in kurtosis(iris$Sepal.Length): 没有"kurtosis"这个函数`

1	`# kurtosis函数默认以3为阈值，所以需要-3使其以0为阈值`

数据处理函数

运算函数

1
2
3

x = 16.139
n = 3
abs(x)

1	`## [1] 16.139`

sqrt(x)

1	`## [1] 4.017337`

1	`round(x,n)# 保留n位小数`

1	`## [1] 16.139`

1	`ceiling(x) # 向上取整`

1	`## [1] 17`

1	`floor(x) # 向下取整`

1	`## [1] 16`

1	`trunc(x) # 向0取整`

1	`## [1] 16`

sin(x)

1	`## [1] -0.4178129`

cos(x)

1	`## [1] -0.9085331`

tan(x)

1	`## [1] 0.4598764`

1	`log(x,n)`

1	`## [1] 2.531593`

1	`exp(x) # e^x`

1	`## [1] 10211244`

字符函数
- 计算x中字符数量nchar：(x)
- 提取向量中的部分字符：substr(x,start,stop)
- 在x中搜索pattern, 并以文本replacement 将其替换：sub(pattern, replacement, x)
- 在split处分割字符向量x中的元素：strsplit(x, split)
- 连接字符串，分隔符为sep：paste(……, sep=” “)
- 大写：toupper(x)
- 小写：tolower(x)
数据处理函数
- 输入向量，返回排序数据：sort(x, na.last=NA, decreasing=FALSE)
- 输入向量，返回排序后数据的index：order(x, na.last=NA, decreasing=FALSE)
  - x[order(x)] == sort(x)
- 返回数值对应的秩次/排位：rank(x, na.last=TRUE, ties.method=c(“average”, “first”, “random”, “max”, “min”))
  - sort(x)[rank(x, ties.method= “first”)] == x
  - ties.method：处理相同值的方法
- 按某一列的值给全表格排序：csv_sorted_age = csv[order(csv$age),]
- 横向合并数据框
  merge(x, y, by, by.x, by.y, all, all.x, all.y)
- 抽样函数
  - sample(x, size, replace=FALSE, prob=NULL) x, 数据； size, 抽取样本数； replace, 放回/不放回； prob, 权重向量
  - seq(from = 1, to = 1, by = ((to - from)/(length.out - 1)), length.out = NULL) from, 开始数字；to, 最大数字； by, 间隔； length.out, 分割成几分
- 数据重塑
  - melt(data, id, measure, variable.name, value.name, na.rm)
    data：要重塑的数据框。
    id.vars：指定那些在重塑后仍然保持为列的变量，通常是标识变量，例如个体的编号、时间等。
    measure.vars：指定那些要被重塑的测量变量列，这些列中的数据会被转换为长格式中的两列（变量名和值）。
    variable.name：指定在重塑后包含原测量变量名称的新列的名字，默认是 “variable”。
    value.name：指定在重塑后包含原测量变量值的新列的名字，默认是 “value”。
  - melt()的逆操作：dcast(data, formula, fun.aggregate = NULL, value.var, …)
    formula：指定如何将数据分组和展开为宽格式。格式一般为 row_variables ~ column_variables，其中 row_variables 是将作为宽格式行的变量，column_variables 是将作为宽格式列的变量。
  - melt()+dcast()的联合操作：recast(data, formula, measure.vars, id.vars, fun.aggregate = NULL, value.var, …)
字符串函数
- 查看长度：str_length(x)
- 字符合并：str_c(x,collapse = “ “)
- 去掉空格和TAB：str_trim(x)
- 字符串值排序：str_sort(sample(LETTERS,10))
- 字符串分割：str_split(x,””)
- 大写：str_to_upper(X)
- 小写：str_to_lower(x)
其他进阶操作
- dplyr包：
  - arrange()排序
    升序：arrange(df, asc(age))
    降序：arrange(df, desc(age))
  - filter()选择行filter(df, age > 30, gender == "male")
  - select()选择列
    - 返回指定列select(df, name, age)
    - 返回指定外列select(carData, -vs, -am)
  - mutate()添加新列
    - 添加列mutate(df, age_plus_10 = age + 10)
    - 去除列mutate(carData, drat=NULL, qsec=NULL)
  - summarize()数据汇总summarise(df, avg_age = mean(age))
  - table()定性数据整理，做频数分布表
  - group_by()分组group_by(df, gender)
  - inner_join() 内连接，返回两个数据框中都有的行inner_join(df1, df2, by = "id")
  - left_join() 左连接，返回 x 数据框中的所有行，以及 y 数据框中与 x 匹配的行left_join(df1, df2, by = "id")
  - right_join() 右连接，返回 y 数据框中的所有行，以及 x 数据框中与 y 匹配的行right_join(df1, df2, by = "id")
  - full_join() 全连接，返回两个数据框中的所有行full_join(df1, df2, by = "id")
- 管道%>%
- attach()函数

1	`attach(iris)`

查看变量类型：class()函数

1
2
3

data = data.frame(c(1,2,3),c(1,2,3))
data = as.matrix(data)
class(data)

1	`## [1] "matrix" "array"`

rbind()和cbind()函数
rbind函数用于按行合并矩阵，即将一个矩阵的每一行添加到另一个矩阵的末尾。合并的前提是两个矩阵的列数必须相同。例如，如果有两个矩阵A和B，A有3行4列，B也有3行4列，使用**rbind(A,B)将得到一个6行4列的新矩阵。如果列数不匹配，R会抛出错误。
cbind函数用于按列合并矩阵，即将一个矩阵的每一列添加到另一个矩阵的末尾。合并的前提是两个矩阵的行数必须相同。例如，如果有两个矩阵A和C，A有3行4列，C有3行5列，使用cbind(A,C)**将得到一个3行9列的新矩阵。如果行数不匹配，R会抛出错误。

a <- matrix(1:12, 3, 4)
b <- matrix(1:12, 3, 4)
y <- rbind(a, b)
print(y)

##      [,1] [,2] [,3] [,4]
## [1,]    1    4    7   10
## [2,]    2    5    8   11
## [3,]    3    6    9   12
## [4,]    1    4    7   10
## [5,]    2    5    8   11
## [6,]    3    6    9   12

样本、统计量与抽样分布

总体：总体(population)是指按研究目的所确定的研究对象中，所有观测个体的某项观测指标取值的集合。总体是服从某种分布的随机变量, 记为 X。
样本：从总体中随机抽取的部分个体，即是样本。每个样本与总体 X 的分布相同。 (抽样的代表性)
- 简单随机样本：随机独立抽取的样本 X1, X2,……,Xn
样本观测值：x1,x2,……,xn

统计量

PPT7.2的4-15页

不含任何未知参数的观测样本的函数

定义：设X1， X2， …， Xn是来自总体 X 的样本,x1， x2， …， xn为其样本值, 则称不含任何未知参数的函数为统计量, 相应实数称为其样本观察(或观测)值
常用统计量
- 样本平均值：$\bar{X} = \frac{1}{n}\sum{X_i}$
- 样本方差（修正）：$S^2 = \frac{1}{n-1}\sum(X_i-\bar{X})^2 = \frac{1}{n-1}(sum{X_i^2-n\bar{X}^2})$
- 样本标准差（修正）：$S = \sqrt{S^2} = \sqrt{\frac{1}{n-1}\sum(X_i-X)^2}$
- 样本标准误：$S_{\bar{X}}=\frac{S}{\sqrt{n}}$
  
  PPT7.2的14-15页：标准差与标准误的辨析
- 样本k阶原点矩：$A_k = \frac{1}{n}\sum{X_i^k},k=1,2,……$
- 样本k阶中心矩：$B_k = \frac{1}{n}\sum(X_i-\bar{X})^k,k=2,3,……$
定理

PPT7.2的10-13页

设总体X的均值为μ，方差为σ^2，X1，……，Xn是总体X的一个样本，则有：$E(\bar{X}) = \mu,D(\bar{X}) = \frac{\sigma^2}{n},E(S^2) = \sigma^2$
若总体X的k阶矩$E(X^k)$记成$\mu_k$存在，则当$n\to\infty$，$A_k\to\mu_k$

抽样分布

PPT7.2的16-21、42-57页

统计量的分布称为抽样分布

定理
- 大样本下，样本均值近似服从正态分布（不要求总体正态）：EXi = EX = μ，DXi = DX = σ^2
- 正态总体的样本均值服从正态分布：$\bar{X} \sim N(\mu,\frac{\sigma^2}{n})$
- 正态总体样本均值与方差的分布：$\bar{X} \sim N(\mu,\frac{\sigma^2}{n})$，$\frac{(n-1)S^2}{\sigma^2} \sim {\chi}^2(n-1)$，$\frac{\bar{X}-\mu}{S/\sqrt{n}} \sim t(n-1)$
- 两独立样本：$\frac{S_1^2/S_2^2}{sigma_1^2/sigma_2^2} \simF(n_1-1,n_2-1)$
  同方差时：$\frac{(\bar{X}-\bar{Y})-(\mu_1-\mu_2)}{S_\omega\sqrt{\frac{1}{n_1}+\frac{2}{n_2}}}$
  其中，$S_\omega^2 = \frac{(n_1-1)S_1^2+(n_2-1)S_2^2}{n_1+n_2-2}$

分布

正态分布

求概率密度

1 2	`a = 72 # a对应的概率密度 dnorm(a, mean=72, sd=15.2)`

1	`## [1] 0.0262462`

区间概率

1 2	`a = 170 # 大于a的概率 pnorm(a, mean=172, sd=5, lower.tail=FALSE)`

1	`## [1] 0.6554217`

1 2	`a = 70 # 大于a的概率 pnorm(a, mean=60, sd=4, lower.tail=FALSE)`

1	`## [1] 0.006209665`

1 2	`a = 170 # 小于a的概率 pnorm(a, mean=172, sd=sqrt(0.08))`

1	`## [1] 7.687299e-13`

1
2
3

a = 4 # 大于a、小于b的概率
b = 10
pnorm(a, mean=1, sd=9, lower.tail=FALSE)-pnorm(10, mean=1, sd=9, lower.tail=FALSE)

1	`## [1] 0.2107861`

已知区间概率，求区间边界

1 2	`a = 0.975 # 概率大于a的区间 qnorm(a, mean=0, sd=1, lower.tail=TRUE)`

1	`## [1] 1.959964`

生成正态总体的样本

1 2	`n = 100 # 样本数 rnorm(n, mean=0, sd=1)`

##   [1]  0.50467868 -0.49163207 -0.43884899  1.10179905  1.99536217  0.64211712
##   [7] -0.08691209  0.72116788 -1.18448775 -0.37059692  0.13703864 -0.31930236
##  [13] -1.01871836  0.01019050 -0.16555060  1.06949034  0.14876524 -0.59551843
##  [19]  0.61369457  0.96441405 -0.39663536  1.09043474  1.23075500 -0.63456115
##  [25] -0.16608714 -0.75793889 -0.32328469 -0.25644695 -0.85971807 -1.83294060
##  [31]  0.53892075 -0.44320198  0.54884601  0.33359645  0.27974270 -0.54318091
##  [37]  0.28101525 -1.46710466 -0.94815968 -1.34588502 -0.85466725 -0.27045963
##  [43]  1.17914320 -0.87250648 -0.69141121 -0.26845721 -1.56518779 -0.11408545
##  [49] -0.24162982  1.10952401 -0.36866151  1.09775125 -0.88393453 -2.14120914
##  [55]  0.48006164 -0.52922710 -0.23827985 -1.45827314 -1.35174748 -0.21767562
##  [61]  0.10449164  0.14103869 -0.97909525 -0.33593721  0.77090022  1.00318046
##  [67]  0.30455796  0.73756327 -0.88370613 -1.00609895  0.98298777  1.23721386
##  [73] -0.68279900 -1.20249314 -0.43900000 -2.30155004  0.65054075  0.56783617
##  [79] -0.82521875  0.78616499  0.96966855  0.08360951  0.08447232  1.57765440
##  [85] -0.27320432  0.44408166 -0.43599298  1.54015523  1.55550025  0.57546912
##  [91] -0.47378658  0.64032989  0.74602902 -0.16675097 -1.94123371  0.50644981
##  [97]  0.90595113 -1.26994551 -1.14191666  0.40864333

t分布

PPT7.2的35页

$\frac{N(0,1)}{\sqrt{\frac{{\chi}^2(n)}{n}}} \sim t(n)$

t分布的密度曲线：关于y轴对称→当$df\to\infty$，为标准正态分布；当df=1，为柯西分布

求概率密度

1
2
3

a = 1 # a对应的概率密度：横坐标→曲线上点
df = 3
dt(a, df)

1	`## [1] 0.2067483`

区间概率

1
2
3

a = 4 # 大于a的概率
df = 3
pt(a, df, lower.tail=FALSE)

1	`## [1] 0.01400423`

1
2
3

a = 4 # 大于a的概率
df = 3
pt(a, df)

1	`## [1] 0.9859958`

a = 4 # 大于a、小于b的概率
b = 10
df = 3
pt(a, df, lower.tail=FALSE)-pt(b, df, lower.tail=FALSE)

1	`## [1] 0.01294003`

已知区间概率，求区间边界/分位点

1
2
3

a = 0.975 # 概率大于a的区间
df = 4
qt(a, df, lower.tail=TRUE)

1	`## [1] 2.776445`

生成正态总体的样本

1
2
3

n = 100 # 样本数
df = 5
rt(n, df)

##   [1]  0.578748722 -1.141025213 -0.256065938 -3.170505942  2.302309922
##   [6] -0.212594614 -0.913242551  0.373528601  0.708768599 -0.072243095
##  [11]  1.499155382  2.022776136  0.009134113  1.600881536  0.715232512
##  [16]  1.063525737  0.553264036  0.143787658  0.057013973  0.249565963
##  [21] -0.640127797 -5.400871115 -0.282034434  0.081255883 -0.242620221
##  [26] -1.706534844  0.460011531 -1.011023681  0.202781803 -1.342483687
##  [31] -0.721754325  0.606157029 -0.607926663 -1.414039512 -0.631158974
##  [36] -0.022346507 -1.790270485  0.667483901  3.951986511  0.367086689
##  [41] -0.030356251 -0.143371501 -1.445657063 -3.855063127  0.174582901
##  [46] -2.419371973  0.506856199  0.849065328  3.008272184 -0.693365559
##  [51]  0.614054352  0.278026556  0.596659912 -0.256291736  5.117292871
##  [56]  0.820699536  0.642524805  0.569704460 -1.830137616  0.457640766
##  [61]  0.183694038 -1.671548525  3.242510765 -0.802718444 -0.323594736
##  [66] -0.254970204 -0.974791967  1.437362739  0.803818866 -0.837508426
##  [71] -0.537972789 -2.095863108  1.429004237  2.916130548  2.011304591
##  [76]  0.051172441  1.008745984 -0.193641295 -2.670538671  3.420001525
##  [81] -1.354037557  1.170810474 -1.137243095  0.030454648  1.808963244
##  [86]  0.136224243  2.119997018  0.250373124  1.180425991 -1.293751225
##  [91] -0.275995313 -3.317088660 -1.081114749 -0.074153765 -0.797602454
##  [96] -0.926499835 -0.116315147  0.105990155  1.738033587 -0.622219714

χ2分布

PPT7.2的22-33页

设X1,X2,……，Xn是来自总体N(0,1)的独立样本，则称统计量$\chi^2 = X_1^2+X_2^2+……+X_n^2$服从自由度为n的$\chi^2$分布，记为$\chi^2(n)$
自由度：指$\chi^2 = X_1^2+X_2^2+……+X_n^2$中包含独立变量的个数
$N^2(0,1) \sim {\chi}^2(1)$

求概率密度

1
2
3

a = 1 # a对应的概率密度：横坐标→曲线上点
df = 3
dchisq(a, df)

1	`## [1] 0.2419707`

区间概率

1
2
3

a = 4 # 大于a的概率
df = 3
pchisq(a, df, lower.tail=FALSE)

1	`## [1] 0.2614641`

1
2
3

a = 4 # 大于a的概率
df = 3
pchisq(a, df)

1	`## [1] 0.7385359`

a = 4 # 大于a、小于b的概率
b = 10
df = 3
pchisq(a, df, lower.tail=FALSE)-pchisq(b, df, lower.tail=FALSE)

1	`## [1] 0.242898`

已知区间概率，求区间边界

1
2
3

a = 0.975 # 概率大于a的区间
df = 4
qchisq(a, df, lower.tail=TRUE)

1	`## [1] 11.14329`

生成正态总体的样本

1
2
3

n = 100 # 样本数
df = 5
rchisq(n, df)

##   [1]  6.2045639  0.8201045  5.0274797  3.1781748  1.1967375 12.7434126
##   [7]  6.2562385  7.5884120  2.7442380  2.9835652 11.4024275  2.5071858
##  [13]  5.5411296  4.2059313  3.9712106  9.5584289 15.5346560  0.9010377
##  [19]  1.7976044  5.1205235  8.9288571  5.9308736  7.5244890  0.6474670
##  [25]  6.4683069  5.0979525 10.1612743  4.7676458  5.4018525  7.2027538
##  [31]  0.9774393  0.9755740  3.9998034  4.3681522  4.3123456 12.8406779
##  [37]  2.3263856  2.2858470  3.8642551  2.3880276  1.2610977  6.1217418
##  [43]  3.6875516  8.0663853 10.0392113  2.6281862  7.9198114  1.6896444
##  [49]  2.8337125  7.7149138  8.2801813  8.7503187  6.2403959  1.1958872
##  [55]  6.6339886  5.9898249  3.3336869  4.5159234  3.7150096  2.2191088
##  [61]  9.6946609  7.3784237 16.1185272  2.3895160  5.6126656  9.7140786
##  [67]  8.2130110  3.2341901  1.1638975  3.6576690  3.1274317  6.5338405
##  [73]  7.5080975  8.8970139  4.3356547  0.8576385  6.7492037  2.8696953
##  [79]  6.4046086  3.5443664  3.2887537  8.1413036  1.9573152  1.4968379
##  [85]  6.5435265  3.5039394  6.7670428  5.7584139  2.5280248  3.2082602
##  [91]  2.6532309  2.7260452  4.0386757 10.5077578  2.5672436  5.4069874
##  [97]  6.7672092  9.4662282  1.1921961  4.5083885

F分布

$\frac{\frac{\chi^2(n_1)}{n_1}}{\frac{\chi^2(n_2)}{n_2}} \sim F(n_1,n_2)$
F分布的密度函数不是对称分布，有关系$F_{1-\alpha}(n,m) = \frac{1}{F_{\alpha}(m,n)}$

求概率密度

a = 1 # a对应的概率密度：横坐标→曲线上点
n1 = 2
n2 = 5
df(a,n1,n2)

1	`## [1] 0.3080008`

区间概率

a = 4 # 大于a的概率
n1 = 2
n2 = 5
pf(a,n1,n2,lower.tail=FALSE)

1	`## [1] 0.09174167`

a = 4 # 大于a的概率
n1 = 2
n2 = 5
pf(a,n1,n2)

1	`## [1] 0.9082583`

a = 4 # 大于a、小于b的概率
b = 10
n1 = 2
n2 = 5
pf(a,n1,n2,lower.tail=FALSE)-pf(b,n1,n2,lower.tail=FALSE)

1	`## [1] 0.07385312`

已知区间概率，求区间边界

a = 0.975 # 概率大于a的区间
n1 = 2
n2 = 5
qf(a,n1,n2,lower.tail=TRUE)

1	`## [1] 8.433621`

生成正态总体的样本

n = 100 # 样本数
n1 = 2
n2 = 5
rf(n, n1,n2)

##   [1]  1.69278371  0.16096309  1.54389441  2.10801557  4.25776558  4.43761066
##   [7] 15.43452835  1.45460478  1.01492928  0.78840775  1.79769939  1.38416368
##  [13]  0.45127317  1.72166643  1.46418133  0.09622615  2.72189072  1.27394302
##  [19]  0.43841567  0.84136556  1.32580550  2.36310646  0.35754291  7.21410013
##  [25]  0.12298536  0.22852830  2.15452588  0.25624828  0.90496477  1.16480374
##  [31]  0.40643084  2.54946501  1.48300702  1.34106811  0.15423003  0.80976918
##  [37]  1.41510422  0.56244561  1.03263665  0.06011787  2.79673698  0.95300817
##  [43]  3.30432145  0.52348463  1.15475573  0.09936875  0.95156828  0.52941999
##  [49]  0.08955415  0.03581757  0.61535978  0.51312807  0.85218589  0.82446437
##  [55]  0.43426724  0.94619987 42.15151770  0.04987677  0.71830400  0.41511948
##  [61]  3.29840181  0.70223562  0.60372711  0.53299506  3.23434781  8.13391978
##  [67]  0.81815374  0.66741564  1.28447298  2.30047025  0.28670049  0.83797160
##  [73]  0.42568775  0.87654478  0.27922448  0.21348101  2.77331242  0.47723637
##  [79]  2.55523798  0.55044106  0.54920483  0.19218782  0.07055777  1.58495316
##  [85]  0.37981103  1.55886596  0.03873870  0.38584984  0.31039453  2.11254536
##  [91]  7.97981917  0.67940267  0.05789905  1.60253601  0.13736147  0.11521427
##  [97]  0.51415643  0.36399651  1.53756473  0.22682322

分布之间的关系

N(0,1)^2 = χ2(1)
t(n)^2 = F(n-1,1)

估计量优劣的标准

PPT12的12-

无偏性/可信度

估计值在未知参数真值附近摆动，它的期望值等于未知参数的真值，无偏估计不唯一
$E(\hat{\theta}-\theta) = 0$

样本均值与方差的无偏估计：见【数据概括性度量】/【统计量】

有效性/精度

当样本量相同时,方差小的无偏估计量更为有效
$D(\hat{\theta_1})<D(\hat{\theta_2})$

均值的无偏估计
样本均值为总体均值μ的无偏估计；在所有加权均值中，样本均值是总体均值μ 的最有效估计

点估计

矩估计/点估计

PPT12的1-8页

定义：用样本矩去替换总体矩，用样本矩的函数去替换相应的总体矩的函数
理论依据：切比雪夫L.L.N.
样本 m 阶矩依概率收敛到理论 m 阶矩
样本矩→总体矩
$A_1 = \frac{1}{n}\sum{X_j} = \sum{x_i\frac{f_i}{n}} \to \mu_1 = EX = \sum{x_ip_i}$
$A_2 = \frac{1}{n}\sum{X_i^2} \to \mu_2 = E(X)^2 = \sum{x_i^2p_i}$
$A_m = \frac{1}{m}\sum{X_m^2} \to \mu_m = E(X)^2 = \sum{x_i^mp_i}$
- 均值：$\mu = EX\to \hat{\mu} = \var{X}$
- 方差：$\sigma^2 = E(X^2)-E^2(X) \to \hat{\sigma}^2 = \frac{1}{n}\sum{(X_i-\bar{X})^2}$

极大似然估计

PPT12的8-11页

思想：样本来自于总体，样本应该很好地反映总体的特征。设总体含有待估计的参数θ，它可以取很多值，应该在θ的一切可能取值之中选出一个使得样本观测值出现的概率为最大的 θ 值来作为估计值

区间估计：置信区间

显著性水平与置信区间
- 置信区间长度 = $2*z_{\alpha/2}(\frac{\sigma}{\sqrt{n}})$
- α = 0.05 对应95%置信区间
- 置信水平是指构造总体参数的多个样本区间中，包含总体参数的区间占总数之比，一般用1-α表示，置信区间越大，置信水平越高
- 单尾和双尾：只影响区间范围，区间线下面积对应显著性水平
  95%的单尾检验：(5%,100%)或(0%,95%)
  95%的双尾检验：(2.5%,97.5%)
  p值的变化：双尾检验的p值是单尾检验的p值的两倍
解释

PPT12的50-51页

参数μ的95%CI的含义或解释：设进行100次重复随机抽样，每次抽样计算参数μ的95%置信区间，最终获得100个μ的置信区间，那么，其中约有95个的区间包含参数μ的真实值μ0

置信区间度量的是在某一置信水平下，样本统计值与
总体参数值间的误差范围

性质
- 置信区间是一个随机区间，每次抽样，区间有所不同，有时置信区间不包含真实的参数
- 置信度 α 固定时，置信区间长度越短，其区间估计的精度越高
- 实际应用中，通常给定一定的置信度，求尽可能短的置信区间
- 置信区间与置信水平和样本量等因素有关
  - 在置信水平固定的情况下，样本量越多，置信区间越窄
  - 在样本量相同的情况下，置信水平越高，置信区间越宽

单个正态总体

求上下界

# N分布
# 单侧
miu = 0
sigma2 = 9
a = 0.05
qnorm(a,miu,sqrt(sigma2),lower.tail = FALSE)# 上界

1	`## [1] 4.934561`

1	`qnorm(a,miu,sqrt(sigma2))# 下界`

1	`## [1] -4.934561`

# 双侧
miu = 0
sigma2 = 1
a = 0.05
qnorm(a/2,miu,sqrt(sigma2),lower.tail = FALSE)

1	`## [1] 1.959964`

1	`qnorm(a/2,miu,sqrt(sigma2))`

1	`## [1] -1.959964`

# t分布
# 单侧
n = 3
a = 0.05
qt(a,n,lower.tail = FALSE)# 上界

1	`## [1] 2.353363`

1	`qt(a,n)# 下界`

1	`## [1] -2.353363`

# 双侧
n = 3
a = 0.05
qt(a/2,n,lower.tail = FALSE)

1	`## [1] 3.182446`

1	`qt(a/2,n)`

1	`## [1] -3.182446`

# 卡方分布
# 单侧
n = 4
a = 0.1
qchisq(a,n,lower.tail = FALSE)# 上界

1	`## [1] 7.77944`

1	`qchisq(1-a,n,lower.tail = FALSE)# 下界`

1	`## [1] 1.063623`

# 双侧
n = 4
a = 0.1
qchisq(a/2,n,lower.tail = FALSE)

1	`## [1] 9.487729`

1	`qchisq(1-a/2,n,lower.tail = FALSE)`

1	`## [1] 0.710723`

# F分布
# 单侧
n = 4
m = 5
a = 0.1
qf(a,n,m,lower.tail = FALSE)# 上界

1	`## [1] 3.520196`

1	`qf(1-a,n,m,lower.tail = FALSE)# 下界`

1	`## [1] 0.2468783`

# 双侧
n = 4
m = 5
a = 0.1
qf(a/2,n,m,lower.tail = FALSE)

1	`## [1] 5.192168`

1	`qf(1-a/2,n,m,lower.tail = FALSE)`

1	`## [1] 0.1598451`

均值

PPT12的25-31页

方差已知

# 给出具体数据
data = c(1,2,3,4)
a = 0.05
sigma2 = 3
n = length(data)
mean = mean(data)
mean - sqrt(sigma2)*qnorm(a/2,lower.tail = FALSE)/sqrt(n)

1	`## [1] 0.8026214`

1	`mean + sqrt(sigma2)*qnorm(a/2,lower.tail = FALSE)/sqrt(n)`

1	`## [1] 4.197379`

# 不给具体数据
a = 0.05
sigma2 = 9
n = 100
mean = 122
mean - sqrt(sigma2)*qnorm(a/2,lower.tail = FALSE)/sqrt(n)

1	`## [1] 121.412`

1	`mean + sqrt(sigma2)*qnorm(a/2,lower.tail = FALSE)/sqrt(n)`

1	`## [1] 122.588`

# 给出区间长度，求样本量n
a = 0.05
length = 1
sigma2 = 9
(2*sqrt(sigma2)*qnorm(a/2,lower.tail = FALSE)/length)^2

1	`## [1] 138.2925`

方差未知

data = c(1,2,3,4)
a = 0.05
n = length(data)
mean = mean(data)
sd = sd(data)
mean - sd*qt(a/2,n-1,lower.tail = FALSE)/sqrt(n)

1	`## [1] 0.4457397`

1	`mean + sd*qt(a/2,n-1,lower.tail = FALSE)/sqrt(n)`

1	`## [1] 4.55426`

方差

PPT12的32页

data = c(5.3,6.6,5.2,3.7,4.9)
a = 0.10
n = length(data)
var = var(data)
(n-1)*var(data)/qchisq(a/2,n-1,lower.tail = FALSE)

1	`## [1] 0.4523738`

1	`(n-1)*var(data)/qchisq(1-a/2,n-1,lower.tail = FALSE)`

1	`## [1] 6.038921`

两个正态总体

均值选择检验思路
- 配对设计→单样本\配对t检验
- 成组设计→
  - 总体方差未知→方差齐性检验F检验→
    - 方差齐性→t检验
    - 方差不齐性→t’检验
  - 总体方差已知→u检验

均值差

方差已知

PPT12的39页

data1 = c(1,2,3,4)
data2 = c(2,3,4,5)
a = 0.05
sigma21 = 3
sigma22 = 4
n = length(data1)
m = length(data2)
mean = mean(data1)-mean(data2)
mean - qnorm(a/2,lower.tail = FALSE)*sqrt(sigma21/n+sigma22/m)

1	`## [1] -3.592789`

1	`mean + qnorm(a/2,lower.tail = FALSE)*sqrt(sigma21/n+sigma22/m)`

1	`## [1] 1.592789`

方差未知
- 方差相同
PPT12的40页

data1 = c(1,2,3,4)
data2 = c(2,3,4,5)
a = 0.05
sigma2 = 3
n = length(data1)
m = length(data2)
mean = mean(data1)-mean(data2)
var1 = var(data1)
var2 = var(data2)
Sw = sqrt(((n-1)*var1+(m-1)*var2)/(n+m-2))
mean - qt(a/2,n+m-2,lower.tail = FALSE)*sqrt(1/n+1/m)*Sw

1	`## [1] -3.233715`

1	`mean + qt(a/2,n+m-2,lower.tail = FALSE)sqrt(1/n+1/m)Sw`

1	`## [1] 1.233715`

方差不同

PPT12的41页

data1 = c(1,2,3,4)
data2 = c(2,3,4,5)
a = 0.05
n = length(data1)
m = length(data2)
mean = mean(data1)-mean(data2)
var1 = var(data1)
var2 = var(data2)
v = (var1/n+var2/m)^2/(var1^2/(n^2*(n-1))+var2^2/(m^2*(m-1)))
mean - qt(a/2,v,lower.tail = FALSE)*sqrt(var1/n+var2/m)

1	`## [1] -3.233715`

1	`mean + qt(a/2,v,lower.tail = FALSE)*sqrt(var1/n+var2/m)`

1	`## [1] 1.233715`

方差比

PPT12的42页

均值已知

data1 = c(1,2,3,4)
data2 = c(2,3,4,5)
a = 0.05
n = length(data1)
m = length(data2)
miu1 = 2
miu2 = 3
var1 = sum((data1 - miu1)^2)/n
var2 = sum((data2 - miu2)^2)/m
var1/(var2*qf(a/2,n,m,lower.tail = FALSE))

1	`## [1] 0.1041175`

1	`var1/(var2*qf(1-a/2,n,m,lower.tail = FALSE))`

1	`## [1] 9.60453`

均值未知

PPT12的43页

data1 = c(1,2,3,4)
data2 = c(2,3,4,5)
a = 0.05
n = length(data1)
m = length(data2)
var1 = var(data1)
var2 = var(data2)
var1/(var2*qf(a/2,n-1,m-1,lower.tail = FALSE))

1	`## [1] 0.06477027`

1	`var1/(var2*qf(1-a/2,n-1,m-1,lower.tail = FALSE))`

1	`## [1] 15.43918`

非正态总体

二项分布

PPT12的44-47页

n>=30
- 单个总体

n = 60 # 尝试次数
m = 40 # 成功次数
a = 0.05
est_p = m/n
est_p - qnorm(a/2,lower.tail = FALSE)*sqrt(est_p*(1-est_p)/n)

1	`## [1] 0.5473871`

1	`est_p + qnorm(a/2,lower.tail = FALSE)sqrt(est_p(1-est_p)/n)`

1	`## [1] 0.7859463`

1 2	`# prop.test()函数 prop.test(x = m, n = n, conf.level = a, correct = FALSE)# correct = TRUE：进行连续性校正`

## 
## 	1-sample proportions test without continuity correction
## 
## data:  m out of n, null probability 0.5
## X-squared = 6.6667, df = 1, p-value = 0.009823
## alternative hypothesis: true p is not equal to 0.5
## 5 percent confidence interval:
##  0.6628396 0.6704718
## sample estimates:
##         p 
## 0.6666667

两个总体

n1 = 50
n2 = 60
m1 = 30
m2 = 40
a = 0.05
est_p1 = m1/n1
est_p2 = m2/n2
est_p = est_p1 - est_p2
est_sigma = sqrt(est_p1*(1-est_p1)/n1+est_p2*(1-est_p2)/n2)
est_p - qnorm(a/2,lower.tail = FALSE)*est_sigma

1	`## [1] -0.2474058`

1	`est_p + qnorm(a/2,lower.tail = FALSE)*est_sigma`

1	`## [1] 0.1140724`

n<30

1
2
3

n = 15# 总实验次数
m = 6# 成功次数
binom.test(m,n,conf.level = 0.95)

## 
## 	Exact binomial test
## 
## data:  m and n
## number of successes = 6, number of trials = 15, p-value = 0.6072
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.1633643 0.6771302
## sample estimates:
## probability of success 
##                    0.4

若要求给定的置信区间长度不超过给定的正数𝛿

PPT12的54页

n = 60 # 尝试次数
m = 40 # 成功次数
a = 0.05
delta = 0.1
est_p = m/n
qnorm(a/2,lower.tail = F)*est_p*(1-est_p)/(delta^2/4)

1	`## [1] 174.219`

泊松分布

PPT12的48页

n>=30

n = 60
lambda = 5
a = 0.05
data = rpois(n, lambda)
est_lambda = sum(data)/n
est_lambda - qnorm(a/2,lower.tail = FALSE)*sqrt(est_lambda/n)

1	`## [1] 4.13575`

1	`est_lambda + qnorm(a/2,lower.tail = FALSE)*sqrt(est_lambda/n)`

1	`## [1] 5.230916`

n<30

1 2	`m = 5# 观察到的事件次数 poisson.test(m,conf.level = 0.95)`

## 
## 	Exact Poisson test
## 
## data:  m time base: 1
## number of events = 5, time base = 1, p-value = 0.00366
## alternative hypothesis: true event rate is not equal to 1
## 95 percent confidence interval:
##   1.623486 11.668332
## sample estimates:
## event rate 
##          5

参数假设检验

基本步骤

PPT13的14-21页

两类错误

I类错误：𝐻0正确，却被拒绝（即目前拒绝了H0）——弃真错误
犯错概率：P = 检验水平α
II类错误：𝐻0错误，却没被拒绝（即目前接受了H0）——存伪错误
- 犯错概率：P = β
- 检验的功效：1-β
两者关系：α↓β↑
减少犯第1类错误的概率，会导致犯第2类错误的概率增加
真实差异|μ-μ0|↑、β↓、(1-β)↑
真实差异越大，发现差异机会越大，功效越大
检验原则：在控制犯第一类错误的前提下，尽可能减小犯第二类的概率
+3个应用
1.计算样本量：设定功效，反推样本量
2.评价假设检验方法的优劣：控制α，看β
3.ROC曲线，判断诊断精确性

单个正态总体

t.test 默认两组数据间相互独立（默认参数 paired = FALSE），执行独立样本的 t 检验；
默认方差不相等（默认参数 var.equal = FALSE）；默认的备择假设是双侧的（默认参数 alternative = ‘two.sided’），即执行双侧检验

均值
- H0：均值等于μ0

data = c(1,2,3,4,5)
h0 = 7
a = 0.05
t.test(x = data, y = NULL, alternative ="two.sided", miu = h0, level = 1-a)

## 
## 	One Sample t-test
## 
## data:  data
## t = 4.2426, df = 4, p-value = 0.01324
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.036757 4.963243
## sample estimates:
## mean of x 
##         3

+ H0：均值大于μ0

data = c(1,2,3,4,5)
h0 = 7
a = 0.05
t.test(x = data, y = NULL, alternative ="greater", miu = h0, level = 1-a)

## 
## 	One Sample t-test
## 
## data:  data
## t = 4.2426, df = 4, p-value = 0.006618
## alternative hypothesis: true mean is greater than 0
## 95 percent confidence interval:
##  1.492557      Inf
## sample estimates:
## mean of x 
##         3

+ H0：均值小于μ0

data = c(1,2,3,4,5)
h0 = 7
a = 0.05
t.test(x = data, y = NULL, alternative ="less", miu = h0, level = 1-a)

## 
## 	One Sample t-test
## 
## data:  data
## t = 4.2426, df = 4, p-value = 0.9934
## alternative hypothesis: true mean is less than 0
## 95 percent confidence interval:
##      -Inf 4.507443
## sample estimates:
## mean of x 
##         3

方差

data = c(1,2,3,4,5)
h0 = 3
a = 0.05
n = length(data)
var = var(data)
sigma = h0
df = n-1
chisq2 = df*var/sigma
q = pchisq(chisq2, df)

泊松分布→正态分布：λ检验

x_bar = 1.5
lambda0 = 1
a = 0.05
z = (x_bar - lambda0) / sqrt(lambda0 / n)
p_value = pnorm(z,0,1,lower.tail = FALSE)
up = qnorm(a/2,0,1,lower.tail = FALSE)
down = qnorm(a/2,0,1)
data.frame(z,p_value,up,down)

1 2	`## z p_value up down ## 1 1.118034 0.1317762 1.959964 -1.959964`

两个正态总体

正态性检验

1
2
3

data1 = c(11.3,15.0,15.0,13.5,12.8,10.0,11.0,12.0,13.0,12.3)
data2 = c(14.0,13.8,14.0,13.5,13.5,12.0,14.7,11.4,13.8,12.0)
shapiro.test(data1)

## 
## 	Shapiro-Wilk normality test
## 
## data:  data1
## W = 0.95998, p-value = 0.7857

1	`shapiro.test(data2)# 检测是否正态`

## 
## 	Shapiro-Wilk normality test
## 
## data:  data2
## W = 0.87357, p-value = 0.11

1	`ks.test(data, "pnorm", mean = mean(data), sd = sd(data))# 可以检测多种分布`

## 
## 	Exact one-sample Kolmogorov-Smirnov test
## 
## data:  data
## D = 0.13646, p-value = 0.9998
## alternative hypothesis: two-sided

方差齐次性检验

data1 = c(3,5,7,9,10)
data2 = c(1,2,4,6,8,4)
h0 = 1
a = 0.05
var.test(x=data1, y=data2, ratio = h0,alternative="two.sided",conf.level=1-a )

## 
## 	F test to compare two variances
## 
## data:  data1 and data2
## F = 1.2487, num df = 4, denom df = 5, p-value = 0.7963
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##   0.1690241 11.6937047
## sample estimates:
## ratio of variances 
##           1.248731

1	`bartlett.test(list(data1,data2)) # <0.05则方差是非齐性`

## 
## 	Bartlett test of homogeneity of variances
## 
## data:  list(data1, data2)
## Bartlett's K-squared = 0.049565, df = 1, p-value = 0.8238

配对均值检验
前提：样本差值di独立同分布，且服从正态分布

data1 = c(0.580,0.509,0.5,0.316,0.337,0.517,0.454,0.512,0.997,0.506)
data2 = c(0.84,0.591,0.674,0.632,0.687,0.978,0.750,0.730,1.2,0.87)
a = 0.05
t.test(x=data1, y=data2,df = length(data1),alternative = "two.sided",paired = TRUE,conf.level = 1-a )

## 
## 	Paired t-test
## 
## data:  data1 and data2
## t = -7.926, df = 9, p-value = 2.384e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.3501458 -0.1946542
## sample estimates:
## mean difference 
##         -0.2724

不配对均值检验

1
2
3

x <- c(4.05, 4.18, 5.93, 4.3, 2.41, 5.6, 6.61, 2.98, 5.93, 4.18, 4.05, 3.62)
y <- c(4.54, 4.63, 3.64, 5.07, 6.44, 5.62, 6.14, 4.81, 6.42, 5.35)
shapiro.test(x) # 先检测是否为正态分布

## 
## 	Shapiro-Wilk normality test
## 
## data:  x
## W = 0.93768, p-value = 0.4687

1	`shapiro.test(y)`

## 
## 	Shapiro-Wilk normality test
## 
## data:  y
## W = 0.95199, p-value = 0.6921

1	`bartlett.test(list(a=x,b=y))`

## 
## 	Bartlett test of homogeneity of variances
## 
## data:  list(a = x, b = y)
## Bartlett's K-squared = 1.0346, df = 1, p-value = 0.3091

1	`t.test(x,y, alternative="less", var.equal=TRUE)`

## 
## 	Two Sample t-test
## 
## data:  x and y
## t = -1.6189, df = 20, p-value = 0.06056
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##        -Inf 0.05093173
## sample estimates:
## mean of x mean of y 
##  4.486667  5.266000

方差齐次

data1 = c(3,5,7,9,10)
data2 = c(1,2,4,6,8,4)
a = 0.05
t.test( x=data1, y=data2,alternative = "two.sided",paired = FALSE,conf.level = 1-a )

## 
## 	Welch Two Sample t-test
## 
## data:  data1 and data2
## t = 1.5925, df = 8.199, p-value = 0.149
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.163860  6.430527
## sample estimates:
## mean of x mean of y 
##  6.800000  4.166667

方差不齐次

data1 = rnorm(10,8,2)
data2 = rnorm(10,5,1)
a = 0.05
t.test( x=data1, y=data2,alternative = "two.sided",paired = FALSE,var.equal = FALSE, conf.level = 1-a )

## 
## 	Welch Two Sample t-test
## 
## data:  data1 and data2
## t = 4.8148, df = 10.635, p-value = 0.0005944
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1.62139 4.37315
## sample estimates:
## mean of x mean of y 
##  8.374857  5.377587

非参数假设检验

符号检验与秩检验：
符号检验：binom.test()
Wilcoxon符号秩检验和秩和检验：wilcox.test()
拟合优度检验
Pearson卡方检验：chisq.test()
Kolmogorov-Smirnov检验：ks.test()
列联表检验：
Pearson卡方检验：chisq.test()
Fisher精确检验：fisher.test()

卡方拟合优度检验

单数据是否符合正态分布

1 2	`data = c(1,2,3,4,5) shapiro.test(data)# 检测是否正态`

## 
## 	Shapiro-Wilk normality test
## 
## data:  data
## W = 0.98676, p-value = 0.9672

1	`ks.test(data, "pnorm", mean = mean(data), sd = sd(data))# 可以检测多种分布`

## 
## 	Exact one-sample Kolmogorov-Smirnov test
## 
## data:  data
## D = 0.13646, p-value = 0.9998
## alternative hypothesis: two-sided

两个数据分布是否一样

# H0：大学生选择五个品牌的概率相同
observed = c(315,108,101,32)
expected = c(09/16,3/16,3/16,1/16)
chisq.test(observed,p = expected)

## 
## 	Chi-squared test for given probabilities
## 
## data:  observed
## X-squared = 0.47002, df = 3, p-value = 0.9254

# 双样本ks检验
mu = 70
sd = 10
set.seed(123)
x = round(rnorm(30,mu,sd))
y = round(runif(30, 0, 100))
  # round()函数用于将数值四舍五入到特定的小数位数
ks.test(x, y)

## 
## 	Exact two-sample Kolmogorov-Smirnov test
## 
## data:  x and y
## D = 0.53333, p-value = 0.0001777
## alternative hypothesis: two-sided

看数据是否符合特定分布

# 均匀分布
x = c(64, 68, 86, 71, 71, 87, 75, 57, 63, 66, 82, 74, 74, 71, 64, 88, 75, 50, 77, 65, 59, 68, 60, 63, 64, 53, 78, 72, 59, 83)

ks.test(x, 'punif', 0, 100) # 成绩是否符合从0-100的均匀分布？

1 2	`## Warning in ks.test.default(x, "punif", 0, 100): ties should not be present for ## the one-sample Kolmogorov-Smirnov test`

## 
## 	Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  x
## D = 0.50333, p-value = 5.006e-07
## alternative hypothesis: two-sided

1	`# 有相同值的时候，ks可能不准`

提示
chisq.test()适合检验分组变量
ks.test()不需要将样本分组，适合于检测检验连续变量，不适合处理离散变量
tie较多的时候，ks.test()的结果不可靠

列联表检验

注意
n ≥ 40，且每个格子的理论频数≥ 5 时，采用卡方检验chisq.test()
对于n ≥ 40且有格子的理论频数< 5 但都 ≥ 1 时，采用 Yates 连续性矫正的卡方检验chisq.test(correct = T)
对于n ≤ 40，或者n > 40且至少有格子的理论频数≤1 时，采用 Fisher 精确独立性检验fisher.test()
独立性：两变量是否相关

1
2
3

M <- as.table(rbind(c(66,16), c(28,18))) 
dimnames(M) <- list(cancer = c("Y", "N"),smoke = c("Y", "N")) 
chisq.test(M, correct = F)

## 
## 	Pearson's Chi-squared test
## 
## data:  M
## X-squared = 5.8143, df = 1, p-value = 0.0159

齐次性：如算发病率有无变化

# 2*3列联表
M <- as.table(rbind(c(120,25), c(50,27),c(40,22))) 
p_effective = (120+50+40)/(120+25+50+27+40+22)
p_ineffective = (25+27+22)/(120+25+50+27+40+22)
chisq.test(M, correct = F)

## 
## 	Pearson's Chi-squared test
## 
## data:  M
## X-squared = 11.951, df = 2, p-value = 0.00254

1
2
3

# 2*2列联表
M <- as.table(rbind(c(81,83), c(19,32))) # 第一行，第二行
chisq.test(M, correct = F)

## 
## 	Pearson's Chi-squared test
## 
## data:  M
## X-squared = 2.3028, df = 1, p-value = 0.1291

2*2列联表

n>=40，且每个E>=5：卡方检验

1 2	`M <- as.table(rbind(c(12,32), c(22,14))) # 第一行，第二行 chisq.test(M, correct = F)`

## 
## 	Pearson's Chi-squared test
## 
## data:  M
## X-squared = 9.2774, df = 1, p-value = 0.00232

n>=40, 1<E<5：yates连续性校正

a = 30# 1行1
b = 45# 1行2
c = 3# 2行1
d = 50# 2行2
chisq2 = (a+b+c+d)*(abs(a*d-b*c)-0.5*n)^2/((a+b)*(c+d)*(a+c)*(b+d))
chisq2

1	`## [1] 19.06814`

1
2
3

# 或者用chisq.test函数
M <- as.table(rbind(c(81,83), c(19,32))) # 第一行，第二行
chisq.test(M, correct = T)

## 
## 	Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  M
## X-squared = 1.8409, df = 1, p-value = 0.1749

n<=40或E<1：fisher精确检验

1 2	`data = matrix(c(1,2,3,4),ncol = 2,nrow = 2) fisher.test(data)`

## 
## 	Fisher's Exact Test for Count Data
## 
## data:  data
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   0.008512238 20.296715040
## sample estimates:
## odds ratio 
##   0.693793

计算理论频数矩阵

# 定义一个函数来计算理论频数矩阵
calculate_expected_frequencies <- function(observed_matrix) {
  # 检查输入矩阵是否为2x2
  if (nrow(observed_matrix) != 2 || ncol(observed_matrix) != 2) {
    stop("输入的矩阵必须是2x2的列联表")
  }
  
  # 计算行总和、列总和和总和
  row_sums <- rowSums(observed_matrix)
  col_sums <- colSums(observed_matrix)
  total_sum <- sum(observed_matrix)
  
  # 计算理论频数矩阵
  expected_matrix <- matrix(0, nrow = 2, ncol = 2)
  for (i in 1:2) {
    for (j in 1:2) {
      expected_matrix[i, j] <- (row_sums[i] * col_sums[j]) / total_sum
    }
  }
  
  return(expected_matrix)
}

# 输入一个2x2的观测频数矩阵
observed_matrix <- matrix(c(13, 7, 4987, 9993), nrow = 2)

# 调用函数计算理论频数矩阵
expected_matrix <- calculate_expected_frequencies(observed_matrix)

# 输出理论频数矩阵
print("观测频数矩阵：")

1	`## [1] "观测频数矩阵："`

1	`print(observed_matrix)`

1
2
3

##      [,1] [,2]
## [1,]   13 4987
## [2,]    7 9993

1	`print("理论频数矩阵：")`

1	`## [1] "理论频数矩阵："`

1	`print(expected_matrix)`

1
2
3

##           [,1]     [,2]
## [1,]  6.666667 4993.333
## [2,] 13.333333 9986.667

1	`chisq.test(observed_matrix, correct = FALSE)`

## 
## 	Pearson's Chi-squared test
## 
## data:  observed_matrix
## X-squared = 9.037, df = 1, p-value = 0.002646

秩检验

符号检验

# H0：A和B一样好，H1：A和B有差异→统计A比B好的数目m
n = 200
m = 72
p = 0.3
a = 0.05
binom.test(m,n,p,alternative="two.sided", conf.level = 1-a)

## 
## 	Exact binomial test
## 
## data:  m and n
## number of successes = 72, number of trials = 200, p-value = 0.07558
## alternative hypothesis: true probability of success is not equal to 0.3
## 95 percent confidence interval:
##  0.2935047 0.4307139
## sample estimates:
## probability of success 
##                   0.36

# H0:A不优于B，H1：A优于B→统计A比B好的数目m
n = 200 # 总次数
m = 72 #成功次数
a = 0.05
p = 0.3
binom.test(m,n, p=0.3,alternative = "greater", conf.level = 1-a) # alternative = "greater"→A大于B

## 
## 	Exact binomial test
## 
## data:  m and n
## number of successes = 72, number of trials = 200, p-value = 0.03963
## alternative hypothesis: true probability of success is greater than 0.3
## 95 percent confidence interval:
##  0.3035147 1.0000000
## sample estimates:
## probability of success 
##                   0.36

配对符号检验

只看A与B的差异的有无，不看差异的具体大小

# H0：两个总体分布相同
# 小样本
data = c(1,0,0,0,0,0,0)
n_pos = sum(data)
n = 7 # 相等组剔除
p = 0.5
p_value = 0
i = 0
while(i<=n_pos){
  p_value = p_value+choose(n,i)*(p)^i*(1-p)^(n-i)
  i = i+1
}
p_value

1	`## [1] 0.0625`

1
2
3

# 大样本示例
a = 0.05
(n_pos-n/2)/sqrt(n/4)#U统计量

1	`## [1] -1.889822`

1	`-qnorm(a/2)`

1	`## [1] 1.959964`

1	`qnorm(a/2)`

1	`## [1] -1.959964`

# 别听上面的了，不如用wilcox做
before <- c(65, 70, 68, 72, 67, 69, 71, 70, 66, 68)
after <- c(68, 69, 67, 75, 65, 70, 72, 70, 67, 68)
wilcox.test(before, after, paired = TRUE)# 这里是配对的关键

1 2	`## Warning in wilcox.test.default(before, after, paired = TRUE): ## 无法精确计算带连结的p值`

1 2	`## Warning in wilcox.test.default(before, after, paired = TRUE): ## 有0时无法计算精确的p值`

## 
## 	Wilcoxon signed rank test with continuity correction
## 
## data:  before and after
## V = 12, p-value = 0.4291
## alternative hypothesis: true location shift is not equal to 0

分位数点的符号检验

单组样本的符号秩检验（有另一组均值/中位数）

wilcox.test 各参数的默认值：
alternative= ‘two.sided’, mu=0, paired= F, exact= NULL, correct= T,conf.level= 0.95

data = c(44.21,45.30,46.39,49.47,51.05,53.16,53.26,54.37,57.16,67.37,71.05,87.37)
# 目测是否为正态
density = density(data)
plot(density)

1	`hist(data)`

1 2	`# 目测不是正态 wilcox.test(data,mu = 45.30, alternative = "greater")`

1 2	`## Warning in wilcox.test.default(data, mu = 45.3, alternative = "greater"): ## 有0时无法计算精确的p值`

## 
## 	Wilcoxon signed rank test with continuity correction
## 
## data:  data
## V = 65, p-value = 0.00255
## alternative hypothesis: true location is greater than 45.3

配对符号秩检验

求秩

1 2	`data = c(3,6,2,4,7,8,1,3,5) rank(data, na.last = TRUE, ties.method = "average")`

1	`## [1] 3.5 7.0 2.0 5.0 8.0 9.0 1.0 3.5 6.0`

# na.last = TRUE：忽略缺失值；ties.method = "average"：用平均数处理相同值

# 标准化
zscores = (data - mean(data))/var(data)

# H0：两个总体分布相同
# 小样本n<=50
data1 = c(0.75,0.59,0.73,0.68,0.58,0.86)
data2 = c(0.72,0.66,0.74,0.84,0.65,0.85)
wilcox.test(data1,data2,paired = TRUE,alternative = "two.sided")

1 2	`## Warning in wilcox.test.default(data1, data2, paired = TRUE, alternative = ## "two.sided"): 无法精确计算带连结的p值`

## 
## 	Wilcoxon signed rank test with continuity correction
## 
## data:  data1 and data2
## V = 4.5, p-value = 0.2463
## alternative hypothesis: true location shift is not equal to 0

1	`# 大样本n>50`

两个独立样本的秩和检验

# H0：健康男性血液中的MYC基因表达值等同于肝癌男性
# H1：健康男性血液中的MYC基因表达值低于肝癌男性

x <- c(23.6,20.9,29.8,11.8,23.0,24.1,17.6,27.6,32.8,27.5)
y <- c(21.9,12.5,7.3,11.7,12.1,11.5,9.5,9.6,10.5,10.2)
wilcox.test(x, y, alternative = "less")

## 
## 	Wilcoxon rank sum exact test
## 
## data:  x and y
## W = 95, p-value = 0.9999
## alternative hypothesis: true location shift is less than 0

# W的计算
z <- c(x,y)
m <-length(x)
W=sum(rank(z)[1:m])-m*(m+1)/2

方差分析

单因素方差分析

假设

原假设H0：μ1=μ2=μ3=……=μn
备择假设H1：存在μi，与其他组的均值不等
拒绝原假设的意义：有足够证据表明至少有一个组的均值与其他组的均值存在显著差异

F值

F值是用于检验不同组之间均值差异显著性的一个统计量。它是组间均方（MSB）与组内均方（MSW）的比值

分子是组间均方（MSB），即不同组之间的变异程度。它反映了不同组的均值之间的差异。
分母是组内均方（MSW），即每组内部数据的变异程度。它反映了随机误差的大小。

如果组间差异显著，则组间均方会远大于组内均方，导致F值较大。相反，如果组间差异较小或不存在，则F值接近1

F值越大，越有可能拒绝原假设，即认为不同组的均值存在显著差异。通常，将计算得到的F值与临界值进行比较，或者通过p值来判断是否拒绝原假设。如果F值大于临界值或p值小于显著性水平（如0.05），则拒绝原假设

具体操作

整体检验必须是正态、独立、齐性样本
- 正态性检验

1 2	`data = c(1,2,3,4,5) shapiro.test(data)# 检测是否正态`

## 
## 	Shapiro-Wilk normality test
## 
## data:  data
## W = 0.98676, p-value = 0.9672

1	`ks.test(data, "pnorm", mean = mean(data), sd = sd(data))# 可以检测多种分布`

## 
## 	Exact one-sample Kolmogorov-Smirnov test
## 
## data:  data
## D = 0.13646, p-value = 0.9998
## alternative hypothesis: two-sided

齐性检验：各组样本量相等，可以不检验方差齐性

method_A <- c(200,220,210,230,190)
method_B <- c(240,250,230,260,220)
method_C <- c(180,190,170,200,195)
data = c(method_A,method_B,method_C)
method_type <- factor(rep(c("A", "B", "C"), each = 5))
df <- data.frame(weight_gain = data, method_type = method_type)

bartlett.test(weight_gain~ method_type, data = df)# 正态数据

## 
## 	Bartlett test of homogeneity of variances
## 
## data:  weight_gain by method_type
## Bartlett's K-squared = 0.33213, df = 2, p-value = 0.847

1	`library(car)`

1	`## 载入需要的程序包：carData`

1	`leveneTest(weight_gain ~ method_type, data = df)# 不需要数据正态`

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  0.2308 0.7974
##       12

方差分析：aov()函数
aov(formula, data=NULL, projections=FALSE, qr=TRUE,
contrasts=NULL, …)
其中，formula为方差分析公式；data为数据框；projection设置是否返回预测结果；qr设置是否返回QR分解结果；contrasts为公式中一些因子的列表。

公式的表示

符号	用法
~	分隔符号，左边为响应（/因）变量，右边为解释（/自）变量eg：y~A+B+C
+	分隔解释变量
:	表示变量的交互项eg：y $\sim$ A+B+A:B
*	表示所有可能交互项eg：y $\sim$ ABC可展开为：y $\sim$ A+B+C+A:B+A:C+B:C+A:B:C
^	表示交互项达到次数eg：y $\sim$ (A+B+C)^2展开为：y $\sim$ A+B+C+A:B+A:C+B:C
.	表示包含除因变量外的所有变量eg：若一个数据框包括变量y,A、B和C，代码y $\sim$ .可展开为y $\sim$ A+B+C

常用表达式
单因素：y $\sim$ A
双因素：y $\sim$ A*B
用已有数据集

1 2	`data(iris) #导入R自带的Iris数据 library("dplyr"); iris %>% group_by(Species) %>% summarize(sw_mean=mean(Sepal.Width), sw_sd = sd(Sepal.Width))`

1 2	`## ## 载入程序包：'dplyr'`

1
2
3

## The following object is masked from 'package:car':
## 
##     recode

1
2
3

## The following objects are masked from 'package:stats':
## 
##     filter, lag

1
2
3

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## # A tibble: 3 × 3
##   Species    sw_mean sw_sd
##   <fct>        <dbl> <dbl>
## 1 setosa        3.43 0.379
## 2 versicolor    2.77 0.314
## 3 virginica     2.97 0.322

1 2	`# 正态性检验 shapiro.test(iris$Sepal.Width[iris$Species=='setosa'])`

## 
## 	Shapiro-Wilk normality test
## 
## data:  iris$Sepal.Width[iris$Species == "setosa"]
## W = 0.97172, p-value = 0.2715

1	`shapiro.test(iris$Sepal.Width[iris$Species=='versicolor'])`

## 
## 	Shapiro-Wilk normality test
## 
## data:  iris$Sepal.Width[iris$Species == "versicolor"]
## W = 0.97413, p-value = 0.338

1	`shapiro.test(iris$Sepal.Width[iris$Species=='virginica'])`

## 
## 	Shapiro-Wilk normality test
## 
## data:  iris$Sepal.Width[iris$Species == "virginica"]
## W = 0.96739, p-value = 0.1809

1 2	`#方差齐性检验 bartlett.test(Sepal.Width~Species,data=iris)`

## 
## 	Bartlett test of homogeneity of variances
## 
## data:  Sepal.Width by Species
## Bartlett's K-squared = 2.0911, df = 2, p-value = 0.3515

1 2	`aov1 <- aov(Sepal.Length~Species, iris) # Species“因素”，Length是结果变量 aov1`

## Call:
##    aov(formula = Sepal.Length ~ Species, data = iris)
## 
## Terms:
##                  Species Residuals
## Sum of Squares  63.21213  38.95620
## Deg. of Freedom        2       147
## 
## Residual standard error: 0.5147894
## Estimated effects may be unbalanced

1	`summary(aov1)`

##              Df Sum Sq Mean Sq F value Pr(>F)    
## Species       2  63.21  31.606   119.3 <2e-16 ***
## Residuals   147  38.96   0.265                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1
2
3

# 提取精确P值
ann<-anova(aov1)
ann[["Pr(>F)"]][1]

1	`## [1] 1.669669e-31`

1 2	`#不同因素水平两两比较可以采用Tukey Honest Significant Differences (TukeyHSD) 统计方法： TukeyHSD(aov1)`

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Sepal.Length ~ Species, data = iris)
## 
## $Species
##                       diff       lwr       upr p adj
## versicolor-setosa    0.930 0.6862273 1.1737727     0
## virginica-setosa     1.582 1.3382273 1.8257727     0
## virginica-versicolor 0.652 0.4082273 0.8957727     0

用自建数据

# 输入数据：按自变量为组别
method_A <- c(78,82,75,80,85)
method_B <- c(85,88,83,90,87)
method_C <- c(70,72,68,75,73)
method_D <- c(88,92,86,90,89)

# 将数据合并成一个向量，并创建对应的因素变量
data <- c(method_A,method_B,method_C,method_D)
methods_type <- factor(rep(c("A", "B", "C","D"), each = 5))

# 进行单因素方差分析
result <- aov(data ~ methods_type)
summary(result)

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## methods_type  3  910.8  303.60   35.61 2.6e-07 ***
## Residuals    16  136.4    8.53                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

计算SSA、SSE、F值和p值，给每组数据

method_A <- c(200,220,210,230,190)
method_B <- c(240,250,230,260,220)
method_C <- c(180,190,170,200,195)

n = 3 # 组数
m = 5 # 每组的样本数
# 有样本具体数值
data = c(method_A,method_B,method_C)
mean_list = c(mean(method_A),mean(method_B),mean(method_C))
mean = mean(data)
var_list = c(var(method_A),var(method_B),var(method_C))

# 只有样本平均数与方差
mean_list = c(33.62,37.83,35.10,34.37)
var_list = c(5.12,4.87,10.98,7.08)
mean = mean(mean_list)
n = 4
m = 10

SSA = m*sum((mean_list-mean)^2)# 每组样本数不一样的话，mi*对应的组方差
SSE = (m-1)*(sum(var_list))# 每组样本数不一样的话，（mi-1）*对应的组方差
N = m*n
F_value = SSA*(N-n)/(SSE*(n-1))
pF = pf(F_value,df1 = n-1,df2 = N-n,lower.tail = FALSE) # 生成的p值是单边检验的
df_A = n-1
df_E = N-n
MSA = SSA/df_A
MSE = SSE/df_E
data.frame(SSA,df_A,MSA,SSE,df_E,MSE,F_value,pF)

1 2	`## SSA df_A MSA SSE df_E MSE F_value pF ## 1 101.086 3 33.69533 252.45 36 7.0125 4.805039 0.006466355`

计算SSA、SSE、F值和p值，给每组均值

n = 3 # 组数
m = 5 # 每组的样本数
mean_A = 1
mean_B = 2
mean_C = 3
s = 10
s_1 = 13
s_2 = 12
s_3 = 11
mean_list = c(mean_A,mean_B,mean_C)
mean = mean(m*(mean_A+mean_B+mean_C))

SSA = m*sum((mean_list-mean)^2)# 每组样本数不一样的话，mi*对应的组方差
SST = (m*n-1)*s
SSE = (m-1)*(s_1+s_2+s_3)# 每组样本数不一样的话，（mi-1）*对应的组方差
F_value = SSA*(m*n - n)/(SSE*(n-1))
pF = pf(F_value,df1 = n - 1,df2 = m*n - n,lower.tail = FALSE) # 生成的p值是单边检验的
df_A = n - 1
df_E = m*n - n
MSA = SSA/df_A
MSE = SSE/df_E
data.frame(SSA,df_A,MSA,SSE,df_E,MSE,F_value,pF)

1 2	`## SSA df_A MSA SSE df_E MSE F_value pF ## 1 11770 2 5885 144 12 12 490.4167 3.117665e-12`

方差不齐/数据不正态时：KW检验→wilcox检验的推广
- kruskal.test()函数

# 创建数据
yield_A <- c(20, 23, 25, 27, 30)
yield_B <- c(18, 20, 19, 22, 23)
yield_C <- c(28, 30, 32, 35, 38)
# 合并数据
yield <- c(yield_A, yield_B, yield_C)
# 创建组别向量
group <- factor(rep(c("A", "B", "C"), each = 5))# each = 5表示每个组别有5个数据
# 进行KW检验
result <- kruskal.test(yield ~ group)
print(result)

## 
## 	Kruskal-Wallis rank sum test
## 
## data:  yield by group
## Kruskal-Wallis chi-squared = 10.441, df = 2, p-value = 0.005405

1 2	`# 绘制箱线图 boxplot(yield ~ group, main = "不同肥料对作物产量的影响", xlab = "肥料组别", ylab = "产量")`

按数学方法计算

n = 3
m = 5
N = n*m
yield_A <- c(20, 23, 25, 27, 30)
yield_B <- c(18, 20, 19, 22, 23)
yield_C <- c(28, 30, 32, 35, 38)
yield <- rank(c(yield_A, yield_B, yield_C))

R_A =  sum(yield[1:5])
R_B = sum(yield[6:10])
R_C = sum(yield[11:15])
list = c(R_A/m,R_B/m,R_C/m)

H = 12/(N*(N+1))*sum(m*(list-(N+1)/2)^2)
p_value = pchisq(H,df = n-1, lower.tail = FALSE)
data.frame(H,p_value)

1 2	`## H p_value ## 1 10.385 0.005558094`

局部分析：两两t检验→矫正p值（克服I类错误膨胀）

# 创建数据
yield_A <- c(20, 23, 25, 27, 30)
yield_B <- c(18, 20, 19, 22, 23)
yield_C <- c(28, 30, 32, 35, 38)
# 合并数据
yield <- c(yield_A, yield_B, yield_C)
# 创建组别向量
group <- factor(rep(c("A", "B", "C"), each = 5))# each = 5表示每个组别有5个数据
# 进行KW检验
result <- kruskal.test(yield ~ group)
print(result)

## 
## 	Kruskal-Wallis rank sum test
## 
## data:  yield by group
## Kruskal-Wallis chi-squared = 10.441, df = 2, p-value = 0.005405

1	`pairwise.t.test(x = yield, g=group, p.adjust.method = "bonferroni", alternative = "two.sided")`

## 
## 	Pairwise comparisons using t tests with pooled SD 
## 
## data:  yield and group 
## 
##   A       B      
## B 0.16030 -      
## C 0.01225 0.00031
## 
## P value adjustment method: bonferroni

双因素方差分析

1 2	`data("ToothGrowth") str(ToothGrowth)`

## 'data.frame':	60 obs. of  3 variables:
##  $ len : num  4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
##  $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
##  $ dose: num  0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...

1 2	`ToothGrowth$dose <-as.factor(ToothGrowth$dose) ToothGrowth %>% group_by(supp,dose) %>% summarize(tlen=mean(len),tsd=sd(len))`

1 2	## `summarise()` has grouped output by 'supp'. You can override using the ## `.groups` argument.

## # A tibble: 6 × 4
## # Groups:   supp [2]
##   supp  dose   tlen   tsd
##   <fct> <fct> <dbl> <dbl>
## 1 OJ    0.5   13.2   4.46
## 2 OJ    1     22.7   3.91
## 3 OJ    2     26.1   2.66
## 4 VC    0.5    7.98  2.75
## 5 VC    1     16.8   2.52
## 6 VC    2     26.1   4.80

1 2	`# 正态性检验 shapiro.test(ToothGrowth$len[ToothGrowth$supp=='VC'])`

## 
## 	Shapiro-Wilk normality test
## 
## data:  ToothGrowth$len[ToothGrowth$supp == "VC"]
## W = 0.96567, p-value = 0.4284

1	`shapiro.test(ToothGrowth$len[ToothGrowth$supp=='OJ'])`

## 
## 	Shapiro-Wilk normality test
## 
## data:  ToothGrowth$len[ToothGrowth$supp == "OJ"]
## W = 0.91784, p-value = 0.02359

1	`shapiro.test(ToothGrowth$len[ToothGrowth$dose==0.5])`

## 
## 	Shapiro-Wilk normality test
## 
## data:  ToothGrowth$len[ToothGrowth$dose == 0.5]
## W = 0.94065, p-value = 0.2466

1	`shapiro.test(ToothGrowth$len[ToothGrowth$dose==1.0])`

## 
## 	Shapiro-Wilk normality test
## 
## data:  ToothGrowth$len[ToothGrowth$dose == 1]
## W = 0.93134, p-value = 0.1639

1	`shapiro.test(ToothGrowth$len[ToothGrowth$dose==2.0])`

## 
## 	Shapiro-Wilk normality test
## 
## data:  ToothGrowth$len[ToothGrowth$dose == 2]
## W = 0.97775, p-value = 0.9019

1 2	`# 方差齐性检验 bartlett.test(len~supp, data=ToothGrowth)`

## 
## 	Bartlett test of homogeneity of variances
## 
## data:  len by supp
## Bartlett's K-squared = 1.4217, df = 1, p-value = 0.2331

1	`bartlett.test(len~dose, data=ToothGrowth)`

## 
## 	Bartlett test of homogeneity of variances
## 
## data:  len by dose
## Bartlett's K-squared = 0.66547, df = 2, p-value = 0.717

1
2
3

# 方差分析
an1<-aov(len~dose+supp,data=ToothGrowth)
summary(an1)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## dose         2 2426.4  1213.2   82.81  < 2e-16 ***
## supp         1  205.4   205.4   14.02 0.000429 ***
## Residuals   56  820.4    14.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1
2
3

# 对该单因素试验进行方差分析：获取P值
ann1<-anova(an1)
ann1[["Pr(>F)"]][1]

1	`## [1] 1.871163e-17`

# 因素的交互作用：两因素是否独立
an2<-aov(len~dose+supp+dose:supp, data=ToothGrowth)
ann2<-anova(an2)
ann2[["Pr(>F)"]][1]# 后面括号里是第1或2因素的p值

1	`## [1] 4.046291e-18`

1 2	`# 因素互作的显著性：用Tukey Honest Significant TukeyHSD(an2,which="dose:supp")`

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = len ~ dose + supp + dose:supp, data = ToothGrowth)
## 
## $`dose:supp`
##                 diff        lwr         upr     p adj
## 1:OJ-0.5:OJ     9.47   4.671876  14.2681238 0.0000046
## 2:OJ-0.5:OJ    12.83   8.031876  17.6281238 0.0000000
## 0.5:VC-0.5:OJ  -5.25 -10.048124  -0.4518762 0.0242521
## 1:VC-0.5:OJ     3.54  -1.258124   8.3381238 0.2640208
## 2:VC-0.5:OJ    12.91   8.111876  17.7081238 0.0000000
## 2:OJ-1:OJ       3.36  -1.438124   8.1581238 0.3187361
## 0.5:VC-1:OJ   -14.72 -19.518124  -9.9218762 0.0000000
## 1:VC-1:OJ      -5.93 -10.728124  -1.1318762 0.0073930
## 2:VC-1:OJ       3.44  -1.358124   8.2381238 0.2936430
## 0.5:VC-2:OJ   -18.08 -22.878124 -13.2818762 0.0000000
## 1:VC-2:OJ      -9.29 -14.088124  -4.4918762 0.0000069
## 2:VC-2:OJ       0.08  -4.718124   4.8781238 1.0000000
## 1:VC-0.5:VC     8.79   3.991876  13.5881238 0.0000210
## 2:VC-0.5:VC    18.16  13.361876  22.9581238 0.0000000
## 2:VC-1:VC       9.37   4.571876  14.1681238 0.0000058

非参数KW检验

# H0：四种方法效果无差别，H1：四种方法效果有差别
mydata_kw<-data.frame(
  ratio=c(1.6,1.7,1.9,1.9,2.0,
          1.5,1.6,1.6,1.7,2.0,
          1.4,1.6,1.6,1.7,1.9,
          1.1,1.1,1.2,1.3,1.6),
  group=c(rep("way1",5),
          rep("way2",5),
          rep("way3",5),
          rep("way4",5))
)
kruskal.test(ratio~group,data=mydata_kw)

## 
## 	Kruskal-Wallis rank sum test
## 
## data:  ratio by group
## Kruskal-Wallis chi-squared = 10.241, df = 3, p-value = 0.01663

理工信（上）开卷代码

写在开头

概率论基础

事件与运算

概率

随机变量

概念及其分布

数字特征

大数定理

中心极限定理

数据

数据概括性度量

数据处理函数

样本、统计量与抽样分布

统计量

抽样分布

分布

正态分布

t分布

χ2分布

F分布

分布之间的关系

估计量优劣的标准

无偏性/可信度

有效性/精度

点估计

矩估计/点估计

极大似然估计

区间估计：置信区间

单个正态总体

求上下界

均值

方差

两个正态总体

均值差

方差比

非正态总体

二项分布

泊松分布

参数假设检验

基本步骤

两类错误

单个正态总体

两个正态总体

非参数假设检验

卡方拟合优度检验

列联表检验

2*2列联表

秩检验

符号检验

配对符号检验

分位数点的符号检验

单组样本的符号秩检验（有另一组均值/中位数）

配对符号秩检验

两个独立样本的秩和检验

方差分析

单因素方差分析

假设

F值

具体操作

双因素方差分析

非参数KW检验

相关与线性回归

相关系数

皮尔逊相关系数ρ

计算ρ

假设检验

成对p值计算

斯皮尔曼spearman秩相关系数

kendall-τ相关系数

线性回归模型

线性回归的假设检验：对于β=0的原假设

决定系数R^2

补充内容

多元线性回归

非线性回归

非参数拟合

局部回归

逻辑回归

cox回归