统计学上分布有很多,在R中基本都有描述。因能力有限,我们就挑选几个常用的、比较重要的简单介绍一下每种分布的定义,公式,以及在R中的展示。
统计分布每一种分布有四个函数:d――density(密度函数),p――分布函数,q――分位数函数,r――随机数函数。比如,正态分布的这四个函数为dnorm,pnorm,qnorm,rnorm。下面我们列出各分布后缀,前面加前缀d、p、q或r就构成函数名:norm:正态,t:t分布,f:F分布,chisq:卡方(包括非中心) unif:均匀,exp:指数,weibull:威布尔,gamma:伽玛,beta:贝塔 lnorm:对数正态,logis:逻辑分布,cauchy:柯西, binom:二项分布,geom:几何分布,hyper:超几何,nbinom:负二项,pois:泊松 signrank:符号秩,wilcox:秩和,tukey:学生化极差
下面先列举各种分布:
rnorm(n, mean=0, sd=1) 高斯(正态)分布
rexp(n, rate=1) 指数分布
rgamma(n, shape, scale=1) γ分布
rpois(n, lambda) Poisson分布
rweibull(n, shape, scale=1) Weibull分布
rcauchy(n, location=0, scale=1) Cauchy分布
rbeta(n, shape1, shape2) β分布
rt(n, df) t分布
rf(n, df1, df2) F分布
rchisq(n, df) χ 2 分布
rbinom(n, size, prob)二项分布
rgeom(n, prob)几何分布
rhyper(nn, m, n, k) 超几何分布
rlogis(n, location=0, scale=1) logistic分布
rlnorm(n, meanlog=0, sdlog=1)对数正态
rnbinom(n, size, prob)负二项分布
runif(n, min=0, max=1)均匀分布
rwilcox(nn, m, n), rsignrank(nn, n) Wilcoxon分布
注意了,上面的分布都有一个规律,就是所有的函数前面都有r开始,所以呢,如果想获得概率密度,就用d替换r
如果想获取累计概率密度,就用p替换r
如果想获取分位数,就用q替换r
二项分布:
即重复n次独立的伯努利试验。在每次试验中只有两种可能的结果,两种结果发生与否互相对立,并且相互独立,与其它各次试验结果无关,事件发生与否的概率在每一次独立试验中都保持不变,则这一系列试验总称为n重伯努利实验,当试验次数为1时,二项分布服从0-1分布。
公式:P(ξ=K)= C(n,k) * p^k * (1-p)^(n-k)
其中,P是成功的概率,n是n次独立重复实验,k是n次实验k次发生的概率
期望:Eξ=np
方差:Dξ=np(1-p)
二项分布在R中展现:
p=.4
K=200
n=10000
x=rbinom(n,k,p)
hist(x)
进行标准化处理:
mean=k*p
var=k*p*(1-p)
z=(x-mean)/sqrt(var)
hist(z)
绘制密度图
mean=k*p
var=k*p*(1-p)
z=(x-mean)/sqrt(var)
hist(z)
正态分布:
正态曲线呈钟型,两头低,中间高,左右对称因其曲线呈钟形,因此人们又经常称之为钟形曲线。
若随机变量X服从一个数学期望为μ、方差为σ^2的正态分布,记为N(μ,σ^2)
当μ = 0,σ = 1时的正态分布是标准正态分布。
正态分布在R中的展现:
x=rnorm(k, mean=mean,sd=sqrt(var))
hist(x)
泊松分布:
是一种统计与概率学里常见到的离散概率分布,由法国数学家西莫恩·德尼·泊松(Siméon-Denis Poisson)在1838年时发表。
泊松分布的概率函数:
泊松分布的参数λ是单位时间(或单位面积)内随机事件的平均发生率。 泊松分布适合于描述单位时间内随机事件发生的次数。
泊松分布在R中的展现:
par(mfrow=c(2,2),mar = c(3,4,1,1))
lambda=.5
x=rpois(k, lambda)
hist(x)
lambda=1
x=rpois(k, lambda)
hist(x)
lambda=5
x=rpois(k, lambda)
hist(x)
lambda=10
x=rpois(k, lambda)
hist(x)
二项分布与泊松分布:
当二项分布的n很大而p很小时,泊松分布可作为二项分布的近似,其中λ为np。通常当n≧10,p≦0.1时,就可以用泊松公式近似得计算。
par(mfrow=c(3,3),mar = c(3,4,1,1))
k=10000
p=c(.5, .05, .005)
n=c(10,100,1000)
for (i in p){
for (j in n){
x=rbinom(k,j,i)
hist(x)
}}
卡方分布:
若n个相互独立的随机变量ξ₁、ξ₂、……、ξn ,均服从标准正态分布(也称独立同分布于标准正态分布),则这n个服从标准正态分布的随机变量的平方和构成一新的随机变量,其分布规律称为卡方分布(chi-square distribution)。
卡方分布是由正态分布构造而成的一个新的分布,当自由度n很大时,
分布近似为正态分布。
卡方分布在R中的展示:
k=10000
par(mfrow=c(2,2),mar = c(3,4,1,1))
x=rchisq(k,2)
d=density(x)
plot(d)
x=rchisq(k,5)
d=density(x)
plot(d)
x=rchisq(k,100)
d=density(x)
plot(d)
x=rchisq(k,1000)
d=density(x)
plot(d)
F分布:
F分布定义为:设X、Y为两个独立的随机变量,X服从自由度为k1的卡方分布,Y服从自由度为k2的卡方分布,这2 个独立的卡方分布被各自的自由度除以后的比率这一统计量的分布。即: F分布是服从第一自由度为k1,第二自由度为k2的分布。
k=10000
par(mfrow=c(2,2),mar = c(3,4,1,1))
x=rf(k,1, 100)
hist(x)
x=rf(k,1, 10000)
hist(x)
x=rf(k,10, 10000)
hist(x)
x=rf(k,10000, 10000)
hist(x)
t分布:
t分布曲线形态与n(确切地说与自由度v)大小有关。与标准正态分布曲线相比,自由度v越小,t分布曲线愈平坦,曲线中间愈低,曲线双侧尾部翘得愈高;自由度v愈大,t分布曲线愈接近正态分布曲线,当自由度v=∞时,t分布曲线为标准正态分布曲线。
k=10000
par(mfrow=c(2,2),mar = c(3,4,1,1))
x=rt(k,2)
hist(x)
x=rt(k,5)
hist(x)
x=rt(k,10)
hist(x)
x=rt(k,100)
hist(x)
几种分布关系图示:
i2mean=function(x,n=10){
k=length(x)
nobs=k/n
xm=matrix(x,nobs,n)
y=rowMeans(xm)
return (y)
}
par(mfrow=c(5,1),mar = c(3,4,1,1))
#Binomia
p=.05
n=100
k=10000
x=i2mean(rbinom(k, n,p))
d=density(x)
plot(d,main="Binomial")
#Poisson
lambda=10
x=i2mean(rpois(k, lambda))
d=density(x)
plot(d,main="Poisson")
#Chi-Square
x=i2mean(rchisq(k,5))
d=density(x)
plot(d,main="Chi-square")
#F
x=i2mean(rf(k,10, 10000))
d=density(x)
plot(d,main="F dist")
#t
x=i2mean(rt(k,5))
d=density(x)
plot(d,main="t dist")
数理统计
基础知识
统计量
mean(x,trim=0,na,rm=FALSE)——均值,trim去掉x两端观测值的便利,默认为0,即包括全部数据,na.rm=TRUE允许数据中有缺失
weighted.mean(x,<weigth>)——加权平均值,weigth表示对应权值
median——中值
quantile(x,probs=seq(<start>,<end>,<diff>))——计算百分位数,是五数总和的扩展,probs设置分位数分位点,用seq(0,1,0.2)设置,表示以样本值*20%为间隔划分数据。
var()——样本方差(n-1)
sd——样本标准差(n-1)
cov——协方差
cor——相关矩阵
fivenum(x,na.rm=TRUE)——五数总括:中位数,下上四分位数,最小值,最大值
数学函数
sum(x,y,z,na.rm=FALSE)——x+y+z,na.rm为TURE可以忽略掉na值数据
sum(x>4)——统计向量x中数值大于4的个数
rep(“LOVE!”,<times>)——重复times次,rep(1:3,c(1,2,3))表示1个1,2个2,3个3组成的序列
sqrt()——开平方函数
2^2 和 **——“^”幂运算
abs()——绝对值函数
'%%'——表示求余
'%/%'——求商(整数)
exp : 2.71828…
expm1 : 当x的绝对值比1小很多的时候,它将能更加正确的计算exp(x)-1
log : 对数函数(自然对数)
log10 : 对数(底为10)函数(常用对数)
log2 : 对数(底为2)函数
因为10>e>1,常用对数比自然对数更接近横坐标轴x
log1p()——log(1+p),用来解决对数变换时自变量p=0的情况。指数和对数的变换得出任何值的0次幂都是1
特性:对数螺旋图。当图像呈指数型增长时,常对等式的两边同时取对数已转换成线性关系。
sin : 正弦函数
cos : 余弦函数
tan : 正切函数
asin : 反正弦函数
acos : 反余弦函数
atan : 反正切函数
sinh : 超越正弦函数
cosh : 超越余弦函数
tanh : 超越正切函数
asinh : 反超越正弦函数
acosh : 反超越余弦函数
atanh : 反超越正切函数
logb : 和log函数一样
log1px : 当x的绝对值比1小很多的时候,它将能更加正确的计算log(1+x)
gamma : Γ函数(伽玛函数)
lgamma : 等同于log(gamma(x))
ceiling : 返回大于或等于所给数字表达式的最小整数
floor : 返回小于或等于所 给数字表达式的最大整数
trunc : 截取整数部分
round : 四舍五入
signif(x,a) : 数据截取函数 x:有效位 a:到a位为止
圆周率用 ‘pi’表示
crossprod(A,B)——A %*% t(B) ,内积
tcrosspeod(A,B)——t(A) %*% B,外积
%*%——内积,a1b1+a2b2+...+anbn=a*b*cos<a,b>,crossprod(x)表示x与x的内积。||x||2,矩阵相乘
%o%——外积,a*b*sin<a,b>(矩阵乘法,叉积),tcrossprod(x,y)表示x与y的外积。*表示矩阵中对应元素的乘积!
向量内积(点乘)和向量外积(叉乘)
正态分布
dnorm(x,mean=0,sd=1,log=FALSE)——正态分布的概率密度函数
pnorm(x,mean=0,sd=1)——返回正态分布的分布函数·
rnorm(n,mean=0.sd=1)——生成n个正态分布随机数构成的向量
qnorm()——下分为点函数
qqnorm(data)——画出qq散点图
qqline(data)——低水平作图,用qq图的散点画线
qq.plot(<x>,main='')——qq图检验变量是否为正态分布
简单分析
summary()——描述统计摘要,和 Hmisc()包的describe()类似,会显示NA值,四分位距是第1个(25%取值小于该值)和第3个四分位数(75%取值小于该值)的差值(50%取值的数值),可以衡量变量与其中心值的偏离程度,值越大则偏离越大。
table(<datafame>$<var>)——统计datafame数据中属性变量var的数值取值频数(NA会自动去掉!),列联表
table(<data_var_1>, <data_var_2>)——比较两个data_var,<data_var_1>为列,<data_var_2>为行,先列后行!
xtabs(formular,data)——列联表
ftable( table())——三维列联表
prop.table()——统计所占百分比例
prop.table(table(<data_var_1>, <data_var_2>),<int>)——比较两个data_var所占百分比,<int>填1位按行百分计算,2为列计算
margin.table( table(),<int> )——计算列联表的边际频数(边际求和),<int>=1为按列变量
addmargin.table(table(),<int> )——计算列联表的边际频数(边际求和)并求和,<int>=1为按列变量
as.formula(<string>)——转换为一个R公式,<string>是一个字符串
循环时的判断语句:
ifelse(<test>, <yes>, <no>)——if,else的变种,test是判断语句,其中的判断变量可以是一个向量!yes是True时的赋值,no是False时的赋值
hist(<data>,prob=T,xlab='横坐标标题',main='标题',ylim=0:1,freq,breaks=seq(0,550,2))——prob=T表示是频率直方图,在直角坐标系中,用横轴每个小区间对应一个组的组距,纵轴表示频率与组距的比值,直方图面积之和为1;prob位FALSE表示频数直方图;ylim设置纵坐标的取值范围;freq为TRUE绘出频率直方图,counts绘出频数直方图,FALSE绘出密度直方图。breaks设置直方图横轴取点间隔,如seq(0,550,2)表示间隔为2,从0到550之间的数值。
density(<data>,na.rm=T)——概率密度函数(核密度估计,非参数估计方法),用已知样本估计其密度,作图为lines(density(data),col="blue")
ecdf(data)——经验分布函数,作图plot(ecdf(data),verticasl=FALSE,do.p=FALSE),verticals为TRUE表示画竖线,默认不画。do.p=FALSE表示不画点处的记号
假设检验
分布函数
shapiro.test(data)——正态W检验方法,当p值大于a为正态分布
ks.test(x,y)——经验分布的K-S检验方法,比较x与y的分布是否相同,y是与x比较的数据向量或者是某种分布的名称,ks.test(x, rnorm(length(x), mean(x), sd(x))),或ks.test(x,"pnorm",mean(x),sd(x))
chisq.test(x,y,p)——Pearson拟合优度X2(卡方)检验,x是各个区间的频数,p是原假设落在小区间的理论概率,默认值表示均匀分布,要检验其它分布,比如正态分布时先构造小区间,并计算各个区间的概率值,方法如下:
brk<-cut(x,br=c(-6,-4,-2,0,2,4,6,8))#切分区间
A<-table(brk)#统计频数
p<-pnorm(c(-4,-2,0,2,4,6,8),mean(x),sd(x))#构造正态分布函数
p<-c(p[1],p[2]-p[1],p[3]-p[2],p[4]-p[3],p[5]-p[4],p[6]-p[5],p[7]-p[6])#计算各个区间概率值
chisq.test(A,p=p)
正态总体的均值方差
t.test(x,y,alternative=c("two.sided","less","greater"),var.equal=FALSE)——单个正态总体均值μ或者两个正态总体均值差μ1-μ2的区间估计;alternative表示备择假设:two.side(默认)是双边检验,less表示H1:μ<μ0,greater表示H1:μ>μ0的单边检验(μ0表示原假设);当var.equal=TRUE时,则是双样本方差相同的情况,默认为不同
var.test(x,y)——双样本方差比的区间估计
独立性检验(原假设H0:X与Y独立)
chisq.test(x,correct=FALSE)——卡方检验,x为矩阵,dim(x)=c(2,2),对于大样本(频数大于5)
fisher.test()——单元频数小于5,列联表为2*2
相关性检验(原假设H0:X与Y相互独立)
cor.test(x,y,method=c("pearson","kendall","spearman"))——相关性检验,观察p-value小于0.05则相关。method选择相关性检验方法
秩
rank()——秩统计量
cor.test()——秩相关检验:Spearman,Kendall
wilcox.test(x,y=NULL,mu,alternative,paired=FALSE,exact=FALSE,correct=FALSE,conf.int=FALSE)——秩显著性检验(一个样本来源于总体的检验,显著性差异的检验),Wilcoxon秩和检验(非成对样本的秩次和检验),mu是待检测参数,比如中值,paired逻辑变量,说明变量x,y是否为成对数据,exact说民是否精确计算P值,correct是逻辑变量,说明是否对p值采用连续性修正,conf.int是逻辑变量,给出相应的置信区间。
uniroot(f,interval=c(1,2))——求一元方程根的函数,f是方程,interval是求解根的区间内,返回值root为解
optimize()或 optimise()——求一维变量函数的极小点
nlm(f,p)——求解无约束问题,求解最小值,f是极小的目标函数,p是所有参数的初值,采用Newton型算法求极小,函数返回值是一个列表,包含极小值、极小点的估计值、极小点处的梯度、Hesse矩阵以及求解所需的迭代次数等。
显著性差异检验(方差分析,原假设:相同,相关性)
mcnemar.test(x,y,correct=FALSE)——相同个体上的两次检验,检验两元数据的两个相关分布的频数比变化的显著性,即原假设是相关分布是相同的。y是又因子构成的对象,当x是矩阵时此值无效。
binom.test(x,n,p,alternative=c("two.sided","less","greater"),conf.level=0.95)——二项分布,符号检验(一个样本来源于总体的检验,显著性差异的检验)
aov(x~f)——计算方差分析表,x是与(因子)f对应因素水平的取值,用summary()函数查看信息
aov(x~A+B+A:B)——双因素方差,其中X~A+B中A和B是不同因素的水平因子(不考虑交互作用),A:B代表交互作用生成的因子
p.adjust()——P值调整函数
pairwise.t.test(x,g,p.adjust.method="holm")——多重t检验,p.adjust.method是P值的调整方法,其方法由p.adjust()给出,默认值按Holm方法(”holm“)调整,若为”none“,表示P值不做任何调整。双因素交互作用时g=A:B
shapiro.test(x)——数据的正态W检验
bartlett.test(x~f,data)——Bartlett检验,方差齐性检验
kruskal.test(x~f,data)——Kruskal-Wallis秩和检验,非参数检验法,不满足正态分布
friedman.test(x,f1,f2,data)——Friedman秩和检验,不满足正态分布和方差齐性,f1是不同水平的因子,f2是试验次数的因子
常用模型
1、回归模型
lm(y~.,<data>)——线性回归模型,“.”代表数据中所有除y列以外的变量,变量可以是名义变量(虚拟变量,k个水平因子,生成k-1个辅助变量(值为0或1))
summary()——给出建模的诊断信息:
1、数据拟合的残差(Residual standard error,RSE),残差应该符合N(0,1)正态的,值越小越好
2、检验多元回归方程系数(变量)的重要性,t检验法,Pr>|t|, Pr值越小该系数越重要(拒绝原假设)
3、多元R方或者调整R2方,标识模型与数据的拟合程度,即模型所能解释的数据变差比例,R方越接近1模型拟合越好,越小,越差。调整R方考虑回归模型中参数的数量,更加严格
4、检验解释变量x与目标变量y之间存在的依赖关系,统计量F,用p-value值,p值越小越好
5、绘图检验plot(<lm>)——绘制线性模型,和qq.plot误差的正态QQ图
6、精简线性模型,向后消元法
线性回归模型基础
lm(formula=x~y,data,subset)——回归分析,x是因变量(响应变量),y是自变量(指示变量),formular=y~x是公式,其中若是有x^2项时,应把公式改写为y~I(x^2),subset为可选择向量,表示观察值的子集。例:lm(Y ~ X1 + X2 + I(X2^2) + X1:X2, data = data)
predict(lm(y~x),new,interval=“prediction”,level=0.95)——预测,new为待预测的输入数据,其类型必须为数据框data.frame,如new<-data.frame(x=7),interval=“prediction”表示同时要给出相应的预测区间
predict(lm(y~x))——直接用用原模型的自变量做预测,生成估计值
筛选模型自变量
lm.new<-update(lm.sol,sqrt(.)~.)——修正原有的回归模型,将响应变量做开方变换
update(<lm>, .~. - x1)——移除变量x1后的模型
coef(lm.new)——提取回归系数
回归诊断
1、正态性(QQ图)
plot(x,which)——回归模型残差图,which=1~4分别代表画普通残差与拟合值的残差图,画正态QQ的残差图,画标准化残差的开方与拟合值的残差图,画Cook统
norm.test()——正态性检验,p-value>0.05为正态
计量的残差图
residuals()和resid()——残差
rstandard()——标准化残差
rstudent()——学生化残差
influence.measures(model)——model是由lm或者glm构成的对象,对回归诊断作总括,返回列表中包括,广义线性模型也可以使用
anova(<lm>)——简单线性模型拟合的方差分析(确定各个变量的作用)
anova(<lm1>,<lm2>)——比较两个模型(检验原假设为不同)
2、误差的独立性——car包提供Duerbin_Watson检验函数
3、线性——car包crPlots()绘制成分残差图(偏残差图)可以看因变量与自变量之间是否呈线性
4、同方差性——car包ncvTest()原假设为误差方差不变,若拒绝原假设,则说明存在异方差性
5、多重共线性——car包中的vif()函数计算VIF方差膨胀因子,一般vif>2存在多重共线性问题
异常点分析(影响分析)
hatvalues()和hat()——帽子矩阵
dffits()——DFFITS准则
cooks.distance()——Cook统计量,值越大越有可能是异常值点
covratio()——COVRATIO准则
kappa(z,exact=FALSE)——多重共线性,计算矩阵的条件数k,若k<100则认为多重共线性的程度很小;100<=k<=1000则认为存在中等程度或较强的多重共线性;若k>1000则认为存在严重的多重共线性。z是自变量矩阵(标准化,中心化的?相关矩阵),exact是逻辑变量,当其为TRUE时计算精准条件数,否则计算近似条件数。用eigen(z)计算特征值和特征向量,最小的特征值对应的特征向量为共线的系数。
step()——逐步回归,观察AIC和残差平方和最小,广义线性模型也可以使用
add1()——前进法
drop()——后退法
stepAIC(sol,direction="backward")——MASS包,可以实现逐步回归(向前、向后、向前向后)
预测
predict(<sol>,<newdataframe>,level=0.95,interval="prediction")——回归预测,sol是模型,newdataframe是待预测数据框,level设置置信度,interval="prediction"表示结果要计算置信区间
glm(formula,family=binomial(link=logit),data=data.frame)——广义线性模型,logit默认为二项分布族的链接函数,formula有两种输入方法,一种方法是输入成功和失败的次数,另一种像线性模型的公式输入方式
predict(glm(),data.frame(x=3.5),type="response")——预测广义线性回归模型,type=“response”表示结果为概率值,否则为预测值y
inv.logit()——预测值y的反logit,boot包的函数
glmnet()——正则化glm函数,glmnet包,执行结果的行数越前正则化越强。其输出结果的意义是:
1)DF是指明非0权重个数,但不包括截距项。可以认为大部分输入特征的权重为0时,这个模型就是稀疏的(sparse)。
2)%Dev就是模型的R2
3)超参数(lambda)是正则化参数。lambda越大,说明越在意模型的复杂度,其惩罚越大,使得模型所有权重趋向于0。
plot(lm(y~x),which=1:4,caption=c(“Residuals vs Fitted”,“Normal Q-Q plot”,“Scale-Location plot”,“Cook's distance plot”))——画回归模型残差图,which为1表示画普通残差与拟合值的残差图,2表示画正态QQ的残差图,3表示画标准化残差的开方与拟合值的残差图,4表示画Cook统计量的残差图;caption是图题的内容。
avova(sol1,sol2,test="Chisq")——比较模型两个模型,广义线性模型可用卡方检验(分类变量),不拒绝原假设说明两个没有显著差异,即用较少自变量模型就可以。
非线性模型
poly(想,degree=1)——计算正交多现实,x是数值向量,degree是正交多项式的阶数,并且degree<length(x)样本个数,例如建立二次正交式回归模型:lm(y~1+poly(x,2))
nls(formula,data,start)——求解非线性最小二乘问题,formula是包括变量和非线性拟合的公式,start是初始点,用列表形式给出
nlm(f,p)——非线性最小二乘,构造最小目标函数,方程移项2为0,f是极小的目标函数,p是所有参数的初值,采用Newton型算法求极小,函数返回值是一个列表,minimum的值便是极小值,estimate是参数的估计值。例如:
fn<-function(p,x,y){
f<-y-p[1]*exp(p[2]*x)
res<-sum(f^2)
}
nlm.sol<-nlm(fn,p=c(3,-0.1),x,y)
2、回归树
rpart( y ~., <data>)——rpart包,回归树,叶结点目标变量的平均值就是树的预测值。生成一棵树,再做修剪(防止过度拟合),内部10折交叉验证
printcp(<rt>)——查看回归树结果,rt是指rpart()函数的运行结果模型,plotcp(<rt>)以图形方式显示回归树的参数信息
参数如下:
cp——当偏差的减少小于某一个给定界限值,默认0.01
minsplit——当结点中的样本数量小于某个给定界限时,默认20
maxdepth——当树的深度大于一个给定的界限值,默认30
prune(<rt>,cp)——自行设置cp值的建树
snip.rpart(<rt>, c(4,7))——修剪,需要修剪的那个地方的是结点号c(4,7),指出输出树对象来需要修剪的树的结点号
snip.rpart(<rt>)——交互修剪,点击结点,右击结束
3、随机森林
randomForest(y ~., <data>)——组合模型,由大量树模型构成,回归任务采用预测结果的平均值。
4、支持向量机
svm(<formula>,<data>,gamma=1/ncol(<data>),<cost>)——e1071包,回归任务,<gamma>=0.01,<cost>=100违反边际所引入的损失?
5、时间序列分析
ts(<data>, frequency=12, start=(2006,1))——把一个向量转化为时间序列对象,<data>向量,frequency表示频率,start表示时间起始点
decompose(<data>,type)——把时间序列分解成长期趋势和周期性变化,<data>是设置了频率(周期长度)的时间序列数据,type="additive"为累加形式:长期趋势+周期性变化+随机变化;"multiplicative"分解为累乘形式:长期趋势*周期性变化*随机变化。默认使用"additive"累加形式。函数返回值sol<-decompose()中,sol$trend是时间序列趋势,seasonal是季节性周期变化,random是随机误差。
stl(<data>,"per")——分解时间序列,返回值sol<-stl()中,sol$time.series[, "seasonal"]读取周期性序列seasonal,sol$time.series[, "trend"]读取长期趋势trend。误差可以使用sol$time.series[, "remainder"]读取。
增长率:
diff(data,lag=1)——差分,上下做差,lag控制变量上下间隔为1
ring.growth[t]=(data[t]-data[t-1])/data[t-1]——同比增长率,描述指标变化趋势
sam.per.grown[t]=(data[t]-data[t-T])/data[t-T]——环比增长率,分析周期性变化,避免周期性变化给数据分析带来的影响,T一般以周为单位
移动平均:
filter(x, filter, method=c("convolution", "recursive"), side=2,...)——线性过滤函数,x待转化的向量数据,method=convolution(卷积方法):使用x内部样本组成线性模型(系数ai由filter参数设置的,side参数设置卷积方法是单边或者双边),recursive(递归方法):使用y内部样本以及当前阶段的x样本组成线性模型(系数ai由filter设置)y递归[t]=x[t]+sum(ai*y[t-i])。side为1(单边卷积)y卷积[t]=a1*x[t]+...+a(k+1)*x[t-k],side为2(双边卷积)y卷积[t]=a1*x[t+m]+...+a(m+1)*x[t]
指数平滑:
sol<-HoltWinters(<data>)——实现二次平滑和三次平滑指数。
sol.forst<-forecast.HoltWinters(sol, h=12)——预测HoltWinters函数产生的模型的新时间序列,h表示频率?预测未来12个月
plot.forecast(sol.forst, include=10)——绘制预测图,include=10表明绘制预测前10个月的数据和未来12个月的预测数据
ARIMA模型
ymd()——lubridate包,将"年-月-日"格式的字符串转换成日期对象,(可以比较前后时间)
自相关性
cov(data.frame(x,y))——协方差矩阵S
cor(data.frame(x,y))——相关系数矩阵R
rnorm(n,<mean>,<sd>)
arima.sim(n=100,list(ar=,ma=))——模拟100个样本的模拟序列
lag.plot(data,lag=k,do.line=FALSE)——绘制原始数据和k阶滞后的散点图
acf(data,lag.max=16,ci.type="ma")——计算并绘制自相关图,0阶自相关系数是rxx,所以恒等于1。ci.type="ma"主要是慨率acf的标准误的问题,以使acf图等准确。
pacf(data,lag.max=16)——偏自相关图,消除Xt-1,...,Xt-k+1的影响后,研究Xt和Xt-k的相关性。
Box.test(data,type="Ljung-Box",lag=16,fitdf=p+q)——自相关性检验,p-value<0.05,标识数据data具有自相关,fitdf为自由度参数p+q
arima(data,order=c(p,d,q))——计算模型参数并建模,TSA包中,order设置AR过程的阶数p,差分过程的d(用于稳定化)和MA过程的阶数q。当p=d=0时,表示只使用MA过程对序列建模。结果sol<-arima()调用predict(sol,n.ahead=5)$pred进行预测,n.ahead参数用于设置预测新阶段的数据量(未来5个月),predict(...)$se标准误差SE,用于计算预测范围(预测范围=预测值+-置信度(alpha)*标准误差SE。
eacf(data)——根据凸显中三角区域顶点的行坐标和列坐标分别确定ARMA的p和q
norm.test()——正态性检验,p-value>0.05为正态
tsdiag(sol)——绘制模型残差的散点图、自相关图和不同阶数下的Box.test体检验p-value值
模型评估
RMSE(lm,< which>)——qpcR包中计算均方根误差,计算子集subset
聚类分析
dist(x,method=”euclidean“)——计算距离
”euclidean“Euclid距离;
”maximum“——Chebyshev距离;
”manhattan“绝对值(马氏)距离;
“canberra”Lance距离;
“minkowski”Minkowski闵式距离;
“binary”定性变量的距离
scale(x, center = TRUE, scale = TRUE)——中心化与标准化,center是中心化,scale是标准化。(全选:减去均值,再除以标准差)
hclust(d,method=“complete”)——系统聚类,d是又dist构成的距离结构,method是系统聚类的方法(默认为最长距离法)
“single”最短距离法“;
”complete“最长距离法;
”median“中间距离法;
”mcquitty“Mcquitty相似法;
”average“类平均法
”centroid“重心法
”ward“离差平法和法
plot(hclist(),hang=0.1)——谱系图,hang表示谱系图中各类所在的位置,hang取负值时,表示谱系图从底部画起。
as.dendrogram(hclust(),hang=-1)——将hclust得到的对象强制转换为谱系图
plot(x,type=c(”rectangle“,”triangle“),horiz=FALSE)——谱系图,x为as.dendrogram返回的对象,type是指是矩形或是三角形,horiz是逻辑变量,当horiz为TRUE时,表示谱系图水平放置。
as.dist()——将普通矩阵转化为聚类分析用的距离结构
plclust(x,hang=0.1)——谱系图,旧版停用,已被plot替换
rect.hclust(x,k,h,border)——在谱系图(plclust())中标注聚类情况,确定聚类个数的函数,x是由hclust生成的对象,k是类个数;h是谱系图中的阈值,要求分成的各类的距离大于h;border是数或向量,标明矩形框的颜色;例如:rec.hclust(hclust(),k=3)
kmeans(x,centers,iter.max,nstart=1,algorithm)——K均值方法,centers是聚类的个数或者是初始类的中心,iter.max为最大迭代次数(默认为10),nstart是随机集合的个数(当centers为聚类的个数时),algorithm为动态聚类算法,例如:km<-kmeans(scale(data),4,nstart=20),返回值中,size表示各类的个数,means表示各类均值,Clustering表示聚类后分类情况?,可以用sort(kmeans()$cluser)对分类情况排序
主成分分析
princomp() 和 prcomp()——主成分分析,结果的标准差显示每一个主成分的贡献率(成分方差占总方差的比例),返回值loadings每一列代表每一个成分的载荷因子
summary(x,loadings=FALSE)——提取主成分的信息,x是princomp()得到的对象,loadings是逻辑变量,为TRUE表示显示主成分分析原始变量的系数,False则不显示。返回表中,Standard deviation是标准差,即方差或lambda的开方,Proportion of Variance表示方差的贡献率,Cumulative Proportion表示累积贡献率。
loadings(x)——显示主成分或因子分析中loadings载荷的内容,主成分是对应割裂,即正交矩阵Q;因子分析中是载荷因子矩阵。x是princomp()或者factanal()得到的对象。
predict(x,newdata)——预测主成分的值,x是由princomp()得到的对象,newdata是由预测值构成的数据框,当newdata为默认值时预测已有数据的主成分值。例如predict(<pca>)[,1]——用主成分的第一列作为原有数据的预测结果
screeplot(x,type=c("barplot",”lines“))——主成分的碎石图,确定主成分维数的选择,x是由princomp()得到的对象,type是描述画出的碎石图的类型,”barplot“是直方图,”lines“是直线图。
biplot(x,choices=1:2,scale=1)——画关于主成分的散点图和原坐标在主成分下的方向,x是由princomp()得到的对象,choices选择主成分,默认为第1、2主成分
factanal(x,factor,covmat=NULL,scores=c("none","regression","Bartlett"),rotation=”varimax“)——因子分析,factors是公因子的个数,covmat是样本协方差和相关矩阵,scores因子得分方法,rotation表示旋转,默认为方差最大旋转
cancor(x,y,xcenter=TRUE,ycenter=TRUE)——典型相关分析,xcenter,ycenter是逻辑变量,为TRUE时做数据中心化