统计计算.docx
- 文档编号:12638725
- 上传时间:2023-04-21
- 格式:DOCX
- 页数:25
- 大小:122.48KB
统计计算.docx
《统计计算.docx》由会员分享,可在线阅读,更多相关《统计计算.docx(25页珍藏版)》请在冰豆网上搜索。
统计计算
一、试用R或Matlab语言对如下矩阵进行正交-三角分解,然后求出其特征根和对应的特征向量。
解:
a<-c(1,1,1,1,-2,-1,2,7,0,2,5,3,-3,1,2,6)
A<-array(a,dim=c(4,4))
qr(A)
eigen(qr(A)$qr)
输出结果:
values
4.333517+0.772598i4.333517-0.772598i-2.333517+1.522577i-2.333517-1.522577i
vectors
[,1][,2][,3]
[1,]0.24662970+0.01321091i0.24662970-0.01321091i-0.92848437+0.00000000i
[2,]-0.88545463+0.00000000i-0.88545463+0.00000000i0.02447355-0.07683699i
[3,]-0.00358256+0.01778396i-0.00358256-0.01778396i-0.13987253+0.31674724i
[4,]0.37414998-0.12104562i0.37414998+0.12104562i0.10542553+0.02015463i
[,4]
[1,]-0.92848437+0.00000000i
[2,]0.02447355+0.07683699i
[3,]-0.13987253-0.31674724i
[4,]0.10542553-0.02015463i
二、试用直接抽样法、复合抽样法、舍选抽样法产生密度为
的随机数,并比较他们的优劣性。
解:
(1)直接抽样法
x=rep(0,100)
for(iin1:
100)
{
r1=runif
(1)
y=r1
x[i]=3*(r1^2)/2
}
x
(2)舍选抽样法
y=seq(-1,1,0.02)
x=rep(0,100)
for(iin1:
100)
{
r=runif
(1)
k=floor(100*r)+1
r1=runif
(1)
a=1.5*(y[k-1])^2
b=1.5*(y[k])^2
d=abs(b-a)/(b+a)
r2=runif
(1)
if(r2<=d&b>a)
x[i]=y[k-1]+(y[k]-y[k-1])*sqrt(r1)
elseif(r2<=d&b x[i]=y[k]-(y[k-1]-y[k])*sqrt(1-r1) else x[i]=y[k-1]+(y[k]-y[k-1])*r1 } x (3)复合抽样法 x=rep(0,100) y=0 for(iin1: 100) { r1=runif (1) r2=runif (1) y=-1+2*r1 while(r2>y^2) { r1=runif (1) r2=runif (1) y=-1+2*r1 } x[i]=y } x 对比: 直接抽样法: 简便直接,但是如果分布函数的反函数不容易求出或者不存在,该方法不适用。 复合抽样法: 若原分布函数能分解成几个常见的分布函数,该方法能有效率的生成随机数。 若不能较为容易分解常见分布函数,该方法会使得求解变得复杂。 舍选抽样法: 简单,应用比直接抽样法更广。 三、用模拟的方法近似计算积分 ,并和精确答案进行比较。 解: n<-100000; y<-rexp(n,rate=1); x<-rexp(n,rate=2); l<-x-y l[l>=0]<-1 l[l<0]<-0 r<-sum(l)/n; r 输出结果: 0.33187 计算结果: = 四、考虑简单的线性模型 1.假定 ,且不同 之间是相互独立的。 取 分别为0,1,2,…,100. (1)取 产生相应的随机数 ,运用最小二乘估计去估计系数 并在坐标系中画出实际数据点和相应的拟合曲线。 解: 程序: e=rnorm(101,0,1) x=seq(0,100) y=10+10*x+e a=rep(1,101) X=matrix(data=c(a,x),ncol=2,nrow=101) B=solve(t(X)%*%X)%*%t(X)%*%y y1=rep(0,101) for(iin1: 101) y1[i]=B[1,1]+B[2,1]*x[i] y1 y plot(y~x) abline(B[1,1],B[2,1]) 结果: 系数估计结果为: [1,]9.942019 [2,]10.003323 (2)取 重复步骤 (1). 解: 程序: e=rnorm(101,0,10) x=seq(0,100) y=10+10*x+e a=rep(1,101) X=matrix(data=c(a,x),ncol=2,nrow=101) B=solve(t(X)%*%X)%*%t(X)%*%y y1=rep(0,101) for(iin1: 101) y1[i]=B[1,1]+B[2,1]*x[i] y1 y plot(y~x) abline(B[1,1],B[2,1]) 结果: 结果系数估计为: [1,]9.872274 [2,]10.012772 2.假定 ,且不同 之间是相互独立的,其中 表示服从位置参数为0,尺度参数为a的Cauchy分布;其他参数含义与 (1)中相同。 对 ,和 两种情况,重复1中类似的操作。 解: (1) 程序: e=rcauchy(101,0,1) x=seq(0,100) y=10+10*x+e a=rep(1,101) X=matrix(data=c(a,x),ncol=2,nrow=101) B=solve(t(X)%*%X)%*%t(X)%*%y y1=rep(0,101) for(iin1: 101) y1[i]=B[1,1]+B[2,1]*x[i] y1 y plot(y~x) abline(B[1,1],B[2,1]) 结果: 系数估计结果为: [1,]6.733459 [2,]10.060441 (2) 程序: e=rcauchy(101,0,100) x=seq(0,100) y=10+10*x+e a=rep(1,101) X=matrix(data=c(a,x),ncol=2,nrow=101) B=solve(t(X)%*%X)%*%t(X)%*%y y1=rep(0,101) for(iin1: 101) y1[i]=B[1,1]+B[2,1]*x[i] y1 y plot(y~x) abline(B[1,1],B[2,1]) 结果: 系数估计结果为: [1,]-427.49891 [2,]19.50897 3.把1、2中的参数 的估计改用最小一乘估计去估计,然后重复做第1和第2小问. 解: 1、 (1) 程序: x=seq(0,100) e=rnorm(101,0,1) y=10+10*x+e M=rep(0,101) M1=rep(0,101) w=rep(0,101) A=x[50] B=y[50] x1=x-A y1=y-B s=sum(abs(x1))/2 for(iin1: 101) { if(x1[i]! =0) w[i]=y1[i]/x1[i] else w[i]=NA } bk=sort(w) p=sort.list(w) for(iin1: 101) { k=p[i] M[i]=abs(x1[k]) } M1[1]=M[1] for(iin2: 101) M1[i]=M1[i-1]+M[i] for(iin1: 101) { m1=M1[i] m11=M1[i+1] if(m1 b=bk[i+1] else l=0 } b a=-A*b+B plot(y~x) abline(a,b) 结果: b=10.00062 a=9.916779 (2) 程序: x=seq(0,100) e=rnorm(101,0,10) y=10+10*x+e M=rep(0,101) M1=rep(0,101) w=rep(0,101) A=x[50] B=y[50] x1=x-A y1=y-B s=sum(abs(x1))/2 for(iin1: 101) { if(x1[i]! =0) w[i]=y1[i]/x1[i] else w[i]=NA } bk=sort(w) p=sort.list(w) for(iin1: 101) { k=p[i] M[i]=abs(x1[k]) } M1[1]=M[1] for(iin2: 101) M1[i]=M1[i-1]+M[i] for(iin1: 101) { m1=M1[i] m11=M1[i+1] if(m1 b=bk[i+1] else l=0 } b a=-A*b+B a plot(y~x) abline(a,b) 结果: b=10.16475 a=-16.60505 2、 (1) 程序: x=seq(0,100) e=rcauchy(101,0,1) y=10+10*x+e M=rep(0,101) M1=rep(0,101) w=rep(0,101) A=x[50] B=y[50] x1=x-A y1=y-B s=sum(abs(x1))/2 for(iin1: 101) { if(x1[i]! =0) w[i]=y1[i]/x1[i] else w[i]=NA } bk=sort(w) p=sort.list(w) for(iin1: 101) { k=p[i] M[i]=abs(x1[k]) } M1[1]=M[1] for(iin2: 101) M1[i]=M1[i-1]+M[i] for(iin1: 101) { m1=M1[i] m11=M1[i+1] if(m1 b=bk[i+1] else l=0 } b a=-A*b+B a plot(y~x) abline(a,b) 结果: b=10.00391 a=10.26043 (2) 程序: x=seq(0,100) e=rcauchy(101,0,100) y=10+10*x+e M=rep(0,101) M1=rep(0,101) w=rep(0,101) A=x[50] B=y[50] x1=x-A y1=y-B s=sum(abs(x1))/2 for(iin1: 101) { if(x1[i]! =0) w[i]=y1[i]/x1[i] else w[i]=NA } bk=sort(w) p=sort.list(w) for(iin1: 101) { k=p[i] M[i]=abs(x1[k]) } M1[1]=M[1] for(iin2: 101) M1[i]=M1[i-1]+M[i] for(iin1: 101) { m1=M1[i] m11=M1[i+1] if(m1 b=bk[i+1] else l=0 } b a=-A*b+B a plot(y~x) abline(a,b) 结果: b=2.973483 a=734.9664 4.比较上述情况,你有什么结论。 结论: 对比可知: 当残差的变异程度较小时,最小二乘法和最小一乘法拟合效果都较好。 五、试用所学过的非线性方法求解最小二乘问题 的最小值。 初值取 程序: k<-c(0,0) fr<-function(x){ x1<-x[1] x2<-x[2] 100*(x2-x1*x1)^2+(1-x1)^2 } kk<-optim(k,fr) 最小值: $value [1]3.729052e-09 六、令 与 是独立的 随机变量,令 。 (1)给出模拟 的粗糙估计量; (2)用对偶变量法进行改进; (3)用条件期望法改进粗糙估计量。 解: 令X、Y为服从于b(5,0.45)、且相互独立的随机变量 (1)粗糙估计量 程序: k=100000 M=rep(0,k) x=rep(0,k) y=rep(0,k) n=5 p=0.45 for(iin1: k) { a=0 b=0 j=0 m=0 repeat { r1=runif (1) j=j+1 if(j>=n)break if(r1<=p) a=a+1 elser1=runif (1) } x[i]=a repeat { r2=runif (1) m=m+1 if(m>=n)break if(r2<=p) b=b+1 elser2=runif (1) } y[i]=b M[i]=exp(a*b) } EM=sum(M)/k 结果: θ的粗糙估计量为: 18032.47 (2)对偶变量法 程序: k=100000 M=rep(0,k) x=rep(0,k) y=rep(0,k) n=5 p=0.45 for(iin1: k) { a=0 b=0 j=0 repeat { r=runif (1) j=j+1 if(j>=n)break else { if(r<=p) a=a+1 elseif((1-r)<=p) b=b+1 elser=runif (1) } } x[i]=a y[i]=b M[i]=exp(a*b) } EM=sum(M)/k 结果: θ的估计量为: 21.90046 七、某维修中心在周末只安排一名员工为顾客提供服务。 新来维修的顾客到达后,若已有顾客正在接受服务,则需要排队等候。 假设来维修的顾客到达过程为Poisson过程,平均4人/小时,维修时间服从指数分布,平均需要6分钟。 试用模拟的方法求该系统的平均队长,以及每位顾客花费在系统中的平均时间和员工平均加班时间。 解: q1<-function(lambda,mu,T) { k<-0;wt<-0;wn<-0;ws<-0; tp<-0;nA<-0;n<-0;t<-0 r<-runif (1);tA<--1/lambda*log(r);tD<-Inf repeat { k<-k+1;wt[k]<-t;wn[k]<-n if(tA { ws[k]<-min(tA,tD)-t if(tA { t<-tA;n<-n+1;nA<-nA+1 r<-runif (1);tA<-t-1/lambda*log(r) if(n==1) r<-runif (1);tD<-t-1/mu*log(r) } else { t<-tD;n<-n-1 if(n==0) tD<-Inf else r<-runif (1);tD<-t-1/mu*log(r) } } else { ws[k]<-if(tD==Inf)0elsetD-t if(n>0) { t<-tD;n<-n-1 if(n>0) r<-runif (1);tD<-t-1/mu*log(r) } else tp<-1 } if(tp==1)break } data.frame(Ls=sum(ws*wn)/t,Ws=sum(ws*wn)/nA, Pwait=sum(ws[wn>=1])/t) } q1(lambda=4,mu=10,T=100) 输出结果为: LsWsPwait 10.68916440.16727370.4185664 八、某市政局为了考察该市的6条交通要道是否有同等流量,分别对这6条交通要道进行了为其一个月的观察,得到每条要道的车辆数(单位: 万辆)分别为157,164,165,182,163,169.试分别用Pearson 检验和模拟的方法估计此数据集的 值,判断此6条交通要道的流量是否相同。 解: 1.用Pearson 检验 X=c(157,164,165,182,163,169) chisq.test(X) 输出结果为: X-squared=2.144,df=5,p-value=0.8289 2.模拟的方法 exam5=function(m) { Q=0;N=rep(0,6); Dat=c(157,164,165,182,163,169) L=sum(Dat) Qt=sum(((Dat-L/6)**2)*6/L) for(iin1: m) { Y=sample(1: 6,L,replace=TRUE) N=c(length(Y[Y==1]),length(Y[Y==2]),length(Y[Y==3]), length(Y[Y==4]),length(Y[Y==5]),length(Y[Y==6])) Q[i]=sum((N-L/6)**2*6/L) } P=length(Q[Q>=Qt])/m } p 输出结果为: 0.8317 (3)模拟的p值为0.8317,不拒绝原假设,即有95%的把握认为六条交通要道的流量相同。 与直接用Pearsonχ2检验的结论相同。 九、设 服从二元正态分布,其均值向量未知,协方差阵为 利用Owen(1990)的经验似然方法给出 的均值的经验似然比渐近置信区间(渐近置信水平 )并画出其草图。 如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。 copyright@ 2008-2022 冰点文档网站版权所有 经营许可证编号:鄂ICP备2022015515号-1s)s)s)s)
冰豆网所有资源均是用户自行上传分享,仅供网友学习交流,未经上传用户书面授权,请勿作他用。