田茂再《多元统计分析 》书中代码汇总.docx
- 文档编号:2426604
- 上传时间:2022-10-29
- 格式:DOCX
- 页数:22
- 大小:20.67KB
田茂再《多元统计分析 》书中代码汇总.docx
《田茂再《多元统计分析 》书中代码汇总.docx》由会员分享,可在线阅读,更多相关《田茂再《多元统计分析 》书中代码汇总.docx(22页珍藏版)》请在冰豆网上搜索。
田茂再《多元统计分析》书中代码汇总
《多元统计分析》代码
第1章多元分布
1.5样本分布和极限定理
例1.8
n=5#n=50
X=matrix(0,n,1000)
for(iin1:
1000)
{X[,i]=rbinom(n,1,1/2)
}
XX=(colMeans(X)-1/2)*5^(1/2)
plot(density(XX),type="l",xlab="1000Random.Sample",
ylab="EstimatedandnormalDensity",main="AsymptoticDistribution,
N=5",col="blue",xlim=c(-2,2))
par(new=TRUE)
dat=rnorm(100000,0,1/2)
plot(density(dat),col="red",main="",axes=FALSE)
1.6厚尾分布
curve(dnorm(x),-5,5,xlab="X",ylab="Y",main="distributioncomparison",
xaxp=c(-5,5,2),ylim=c(0,0.4),col="blue")
curve(dcauchy(x,0,1),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,0.4),col="red")
legend("topright",c("Gaussian","Cauchy"),pch=c(16,16),col=c("blue","red")
text.col=c("blue","red"),bty="n")
abline(v=c(-2,-1,1,2))
text(-2.1,0.39,"-2f",col="blue",adj=c(0,-0.1))
text(-1.1,0.39,"-f",col="blue",adj=c(0,-0.1))
text(0.9,0.39,"f",col="blue",adj=c(0,-0.1))
text(1.9,0.39,"2f",col="blue",adj=c(0,-0.1))
1.6.1广义双曲分布
library(Runuran)
#为广义双曲分布产生分布函数目标
distr1<-udghyp(lambda=0.5,alpha=1,beta=0,delta=1,mu=0)
#制造生成器目标;利用PINV(逆)方法
gen1<-pinvd.new(distr1)
#产生大小为10000的样本
x1<-ur(gen1,10000)
#概率密度函数
dx1<-ud(gen1,x1)
#累积分布函数
fx1<-up(gen1,x1)
o1=order(x1)
#为HYD产生分布目标
distr2<-udghyp(lambda=1,alpha=1,beta=0,delta=1,mu=0)
gen2<-pinvd.new(distr2)
x2<-ur(gen2,10000)
dx2<-ud(gen2,x2)
fx2<-up(gen2,x2)
o2=order(x2)
#画图
op<-par(mfrow=c(1,2))
plot(sort(x1),dx1[o1],type="l",xlab="X",ylab="Y",
main="pdfofGH,HYPandNIC",xlim=c(-5,5),
xaxp=c(-5,5,2),ylim=c(0,0.52),col="black")
lines(sort(x2),dx2[o2],type="l",col="red")
lines(sort(x3),dx2[o3],type="l",col="blue")
legend("topright",c("GH","NIC","HYP"),pch=c(16,16,16),
col=c("black","red","blue"),text.col=c("black","red",
"blue"),bty="n")
plot(sort(x1),fx1[o1],type="l",xlab="X",ylab="Y",
main="cdfofGH,HYPandNIC",xlim=c(-5,5),xaxp=
c(-5,5,2),ylim=c(0,1),col="black")
lines(sort(x2),fx2[o2],type="l",col="red")
lines(sort(x3),fx2[o3],type="l",col="blue")
legend("topleft",c("GH","NIG","HYP"),pch=c(16,16,16),col=
c("black","red","blue"),text.col=c("black","red","blue"),
bty="n")
par(op)
1.6.2学生t分布
图1-5
op<-par(mfrow=c(1,2))
curve(dt(x,3),-5,5,xlab="X",ylab="Y",main="pdfoft-distribution",
xaxp=c(-5,5,2),ylim=c(0,0.42),col="black")
curve(dt(x,6),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,0.42),col="blue")
curve(dt(x,30),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,0.42),col="red")
legend("topright",c("t3","t6","t30"),pch=c(16,16,16),col=c("black","blue","red")
text.col=c("black","blue","red"),bty="n")
curve(pt(x,3),-5,5,xlab="X",ylab="Y",main="cdfoft-distribution",
xaxp=c(-5,5,2),ylim=c(0,1),col="black")
curve(pt(x,6),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,1),col="blue")
curve(pt(x,30),-5,5,add=TRUE,xaxp=c(-5,5,2),ylim=c(0,1),col="red")
legend("topleft",c("t3","t6","t30"),pch=c(16,16,16),col=c("black","blue","red")
text.col=c("black","blue","red"),bty="n")
par(op)
图1-6
op<-par(mfrow=c(1,2))
curve(dt(x,1),2.6,3.8,xlab="X",ylab="Y",main="tailcomparison-t-distribution",
xaxp=c(3,3.5,1),ylim=c(0,0.04),col="black")
curve(dt(x,3),2.6,3.8,add=TRUE,xaxp=c(3,3.5,1),ylim=c(0,0.04),col="blue")
curve(dt(x,9),2.6,3.8,add=TRUE,xaxp=c(3,3.5,1),ylim=c(0,0.04),col="red")
curve(dt(x,45),2.6,3.8,add=TRUE,xaxp=c(3,3.5,1),ylim=c(0,0.04),col="deeppink")
curve(dnorm(x),2.6,3.8,add=TRUE,xaxp=c(3,3.5,1),ylim=c(0,0.04),col="grey")
points(x=3.2,y=dt(3.2,1),cex=0.8,pch=1,col="black")
text(x=3.21,y=dt(3.2,1),"t1",cex=0.8,col="black",adj=c(0,-0.1))
points(x=3.2,y=dt(3.2,3),cex=1,pch=1,col="blue")
text(x=3.21,y=dt(3.2,3),"t3",cex=0.8,col="blue",adj=c(0,-0.1))
points(x=3.2,y=dt(3.2,9),cex=1,pch=1,col="red")
text(x=3.21,y=dt(3.2,9),"t9",cex=0.8,col="red",adj=c(0,-0.1))
points(x=3.2,y=dt(3.2,45),cex=1,pch=1,col="deeppink")
text(x=3.21,y=dt(3.2,45),"t45",cex=0.8,col="deeppink",adj=c(0,-0.1))
points(x=3.2,y=dnorm(3.2),cex=1,pch=1,col="grey")
text(x=3.21,y=0,"Gaussian",cex=0.8,col="grey",adj=c(0,-0.1))
curve((abs(x))^(-2),1.2,1.5,xlab="X",ylab="Y",main="tailcomparison-approximation",
xaxp=c(1.2,1.5,6),ylim=c(0,0.7),col="black")
curve((abs(x))^(-4),1.2,1.5,add=TRUE,xaxp=c(1.2,1.5,6),ylim=c(0,0.7),col="black")
curve((abs(x))^(-10),1.2,1.5,add=TRUE,xaxp=c(1.2,1.5,6),ylim=c(0,0.7),col="red")
points(x=1.35,y=(abs(1.36))^(-2),cex=0.8,pch=1,col="black")
text(x=1.36,y=(abs(1.35))^(-2),"t1",cex=0.8,col="black",adj=c(0,-0.1))
points(x=1.35,y=(abs(1.36))^(-4),cex=1,pch=1,col="blue")
text(x=1.36,y=(abs(1.35))^(-2),"t3",cex=0.8,col="blue",adj=c(0,-0.1))
points(x=1.35,y=(abs(1.36))^(-10),cex=1,pch=1,col="red")
text(x=1.36,y=(abs(1.35))^(-10),"t9",cex=0.8,col="red",adj=c(0,-0.1))
par(op)
1.6.3拉普拉斯分布
library(Runuran)
#为拉普拉斯分布产生分布函数目标
distr1<-udlaplace(location=0,scale=1)
#制造生成器目标;利用PINV(逆)方法
gen1<-pinvd.new(distr1)
#产生大小为10000的样本
x1<-ur(gen1,10000)
#密度
dx1<-ud(gen1,x1)
#累积分布函数
fx1<-up(gen1,x1)
o1=order(x1)
distr2<-udlaplace(location=0,scale=1.5)
gen2<-pinvd.new(distr2)
x2<-ur(gen2,10000)
dx2<-ud(gen2,x
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 多元统计分析 田茂再多元统计分析 书中代码汇总 田茂再 多元 统计分析 代码 汇总