统计计算程序重点
- 格式：doc
- 大小：51.50 KB
- 文档页数：5
bootstrap一、function Qb =example_bootstrap(n,sig,B)X = normrnd(15*ones(n,1),sig);%生成均值为15，标准差为sig，服从正态分布的n行1列%随机变量hsig2 = X'*(eye(n) -ones(n,1)*ones(1,n)/n)*X/n;%对上述X的sig2估计值Qb = zeros(B,1);for b = 1:B %重复有放回抽样过程n次（以上述生成的X作为总体）id = unidrnd(n,n,1);%最大值为n 的n行1列随机数。
有放回抽样过程，unidrnd(N)表示产生从1到N所指定的%最大数之间的离散均匀随机整数Xb = X(id);%新样本hsig2b = Xb'*(eye(n) -ones(n,1)*ones(1,n)/n)*Xb/n;%新生成样本的sig2估计值Qb(b) = n*hsig2b/hsig2; %服从自由度为（n-1)的卡方分布endhq = sort(Qb);%排序p0 = ((1:B) - 0.5)'/B;%从小到大的B 个概率值q0 = chi2inv(p0,n-1);%根据概率反查临界值plot(q0,hq,'o');axis([q0(1) q0(B) hq(1) hq(B)])%坐标轴的范围 axis([xmin xmax ymin ymax])二、function [qb,q0] =example_bootstrap_1(n,sig,B)X = normrnd(15*ones(n,1),sig);hsig2 = X'*(eye(n) -ones(n,1)*ones(1,n)/n)*X/n;Qb = zeros(B,1);for b = 1:Bid = unidrnd(n,n,1);Xb = X(id);hsig2b = Xb'*(eye(n) -ones(n,1)*ones(1,n)/n)*Xb/n;Qb(b) = n*hsig2b/hsig2;endqb =quantile(Qb,[0.25;.5;.75]);%y=qu antile(x,50) the median of x%y=quantile(x,[.25 .50 .75]) the quartiles of xq0 = chi2inv([0.25;.5;.75],n-1);第七章交叉验证function [APE] =example_cv(y,x1,x2)DataFull = [y x1 x2];n = length(y);for k = 1:n%%%%%%%% DataPartitioning %%%%%%%%%%%%%%%DataTest = DataFull(k,:);%测试集yTest = DataTest(:,1);%测试集的yXTest = [1DataTest(:,2:3)];%测试集的常数项，x1,x2DataTrain = DataFull;%训练集DataTrain(k,:) = [];%去掉训练集的第k项得到新的训练集yTrain = DataTrain(:,1);%训练集的y,(n-1)行1列XTrain = [ones(n-1,1) DataTrain(:,2:3)];%ones（n-1,1)表示常数项，(n-1)行%%%%%%%%% Define Candidate Models %%%%%%%%%Define Model 1X1 = XTrain(:,1:2);hb1 =regress(yTrain,X1);% OLS estimator 最小二乘法回归hb1 = [hb1;0];%hb1也可写作hb(1);hb(2),是回归得到的常数项和x1的系数py1 = XTest*hb1;% preditionLoss1(k) = (yTest -py1)^2;%prediction error（y-y1bar)^2%Define Model 2X2 =[XTrain(:,1),XTrain(:,3)];hb2 = regress(yTrain,X2); hb2 = [hb2(1);0;hb2(2)]; py2 = XTest*hb2;Loss2(k) = (yTest - py2)^2;%Define Model 3X3 = XTrain;hb3 = regress(yTrain,X3); py3 = XTest*hb3;Loss3(k) = (yTest - py3)^2; endLoss = [Loss1;Loss2;Loss3];APE = mean(Loss,2);%mean(A)和mean(A,1)表示求矩阵A各列的均值；mean(A,2)表示%矩阵A各行的均值更改function [APE] =example_cv01(y,x1,x2)DataFull = [y x1 x2.^2 x1.*x2];n = length(y);for k = 1:n%%%%%%%% DataPartitioning %%%%%%%%%%%%%%%DataTest = DataFull(k,:);yTest = DataTest(:,1);XTest = [1 DataTest(:,2:4)]; DataTrain = DataFull;DataTrain(k,:) = [];yTrain = DataTrain(:,1);XTrain = [ones(n-1,1) DataTrain(:,2:4)];%%%%%%%%% Define Candidate Models %%%%%%%%%Define Model 1X1 = XTrain(:,1:2);hb1 =regress(yTrain,X1);% OLS estimator hb1 = [hb1;0;0];py1 = XTest*hb1;% preditionLoss1(k) = (yTest -py1)^2;%prediction error%Define Model 2X2 =[XTrain(:,1),XTrain(:,3),XTrain( :,4)];hb2 = regress(yTrain,X2); hb2 =[hb2(1);0;hb2(2);hb2(3)];py2 = XTest*hb2;Loss2(k) = (yTest - py2)^2;%Define Model 3X3 =[XTrain(:,1),XTrain(:,2),XTrain( :,3)];hb3 = regress(yTrain,X3); py3 = XTest*[hb3;0];Loss3(k) = (yTest - py3)^2; endLoss = [Loss1;Loss2;Loss3];APE = mean(Loss,2);Jackknifefunction[hsig2,hbias,hsig2BC, hv] = example_JK(y)hsig2 = mean((y - mean(y)).^2);%有偏估计 1/nhv = var(y);%无偏估计 1/(n-1)n = length(y);hsig2jk = zeros(n,1);for k =1:nyjk = y;yjk(k) = [];hsig2jk(k) = mean((yjk -mean(yjk)).^2);endhbias = (n - 1)*(mean(hsig2jk) - hsig2);%biashsig2BC = hsig2 - hbias;%修正后所得值POWER一、样本容量变动function example_6_3(sig)mu0=45:0.05:48;mu0=mu0';nn=[50;100;200];c=norminv(0.95);m=length(nn);pp=zeros(length(mu0),m);for k=1:mn=nn(k);pp(:,k)=1-normcdf(c*sig/sqrt(n), mu0-45,sig/sqrt(n));%表示服从N（mu0-45,%sig/sqrt(n)方，落在点c*sig/sqrt(n)右侧的概率normcdf(x,mu,sigma)end plot(mu0,pp(:,1))grid minorhold onplot(mu0,pp(:,2),'r')plot(mu0,pp(:,3),'g')hold off二、标准差变动function example_6_4(n)mu0=45:0.05:48;mu0=mu0';sig00=[0.5;1;2];c=norminv(0.95);m=length(sig00);pp=zeros(length(mu0),m);for k=1:msig=sig00(k);pp(:,k)=1-normcdf(c*sig/sqrt(n), mu0-45,sig/sqrt(n));endplot(mu0,pp(:,1))grid minorhold onplot(mu0,pp(:,2),'r')plot(mu0,pp(:,3),'g')hold off偏度function [bias1, bias2] =examp_chap3_unbias(n,p,M)s0 = zeros(n,1);s1 = zeros(n,1);%重复M次for k=1:Mx = binornd(1,p,n,1);%服从于b(1,p)的n行1列随机变量s0(k) = x'*(eye(n) -ones(n,n)/n)*x/n;s1(k) = x'*(eye(n) -ones(n,n)/n)*x/(n - 1);%所抽取的样本方差endbias1 = mean(s0) - p*(1-p);%p*(1-p)即为服从b(1,p)的方差bias2 = mean(s1) - p*(1-p);柱状图mu=mean(forearm);v=std(forearm);xp=linspace(min(forearm),max(for earm));yp=normpdf(xp,mu,v);[n,x]=hist(forearm,20);h=x(2)-x(1);bar(x,n/(93*h),1);hold onplot(xp,yp,'r');hold off逆变换法生成随机变量理解用：function[m1,m2,m3,s]=example_4_2(ProbMas s,n)F=cumsum(ProbMass);%thecumulative distribution functionx=zeros(n,1);U=rand(n,1);%生成均匀随机变量for k=1:nx(k)=sum(U(k)>F);%U(k)>F是对于每一个U（k)均与F中的每一项（0.3,0.5,1）%作比较，若为真结果为1，为假结果为0，最后求和并赋值给x(k)endunique(x) 检验m1=mean(x==0);m2=mean(x==1);m3=mean(x==2);s=sum(x==0);通用functionx=example_4_1(ProbMass,n)A=[-1 2 4 5];F=cumsum(ProbMass);x=zeros(n,1);U=rand(n,1);for k=1:nx(k)=A(sum(U(k)>F)+1);%将A(1),A（2），A(3),A(4）赋值给x(k）endpdf:密度 fcdf::累积分布 Fbetapdf(x,t,1/λ)Weibcdf(x,），βαβ-quantile(x,[.25,.5,.75])norminv(0.95) 求临界值tinv(概率，自由度）chi2inv(概率，自由度）Poissinv（概率，lambda)normrnd(mu,sig0,m,n)poissonrnd(lambda,m,n)trnd（v,m,n)chi2rnd(v,m,n)gamma(t)=(t-1)!计算统计和传统统计计算统计是一个有趣的和相对较新的领域内统计。
