财大统计计算复习

load steam                   %加载数据 

n = length(y);

X = [ones(n,1), x', (x').^2]; %X'为列向量

PE1 = zeros(n,1);

PE2 = zeros(n,1);

for k = 1:n

    Ytest = y(k);

    Ytrain = y';

    Ytrain(k) = [];         %表示删掉第y(k)

    Xtest = X(k,:);

    Xtrain = X;

    Xtrain(k,:) = [];

    % evaluate the prediction error of model 1

    X1  = Xtrain(:,1:2);

    hb1 = (X1'*X1)\\(X1'*Ytrain);% evaluate the OLS estimate of beta

    PE1(k) = (Ytest - Xtest*[hb1;0])^2;

    % evaluate the prediction error of model 2

    X2  = Xtrain;

    hb2 = (X2'*X2)\\(X2'*Ytrain);

    PE2(k) = (Ytest - Xtest*hb2)^2;

end

PE = [PE1 PE2];

boxplot(PE)

mean(PE)

q=quantile(PE,3)

第6章 纠偏

1、主程序：

function [s2BC, biasB] = Example610a(X,B)

s2 = var(X,1);                            %var(X,1)表示/n 理论方差, var(X)表示/n-1

n  = length(X);                           %原始数据X的维数

sb2 = zeros(B,1);                         %Bootstrap的样本量为B

for b = 1:B

    id = unidrnd(n,n,1);                  %生成1—n上的n*1列  抽取X的脚标

    Xb = X(id);                           %将对应脚标的X值提取出来赋值给Xb

    sb2(b) = var(Xb,1);                   %用Bootstrap方法抽取的样本的方差

end

biasB = mean(sb2) - s2                    %偏差

s2BC = s2 - biasB                         %纠偏之后的方差

子程序：

clear

n = 200;

sig2 = 1.5;

X = normrnd(2*ones(n,1),sqrt(sig2));

s2 = var(X,1);

B  = 1e4;

[s2BC, biasB] = Example610a(X,B)

2、主程序

function sb2 = Example610b(X,B)

s2 = var(X,1);

n  = length(X);

sb2 = zeros(B,1);

for b = 1:B

    id = unidrnd(n,n,1);

    Xb = X(id);

    sb2(b) = var(Xb,1);

end

子程序：

clear

n = 200;

sig2 = 1.5;

X = normrnd(2*ones(n,1),sqrt(sig2));

s2 = var(X,1);

B  = 1e4;

sb2 = Example610b(X,B);

功效图

%画功效图1-bata

function Pow = Example61a(n)

n=100;

Mu = 119:.01:121;             %不同的均值

M  = 1e4;                     %产生随机向量的个数

Pow = zeros(length(Mu),1);    %功效向量储存向量

for s = 1:length(Mu)          %length表示长度，包含的个数

    mu = Mu(s);               %第Mu的第s个数

    Flag = zeros(M,1);        %储存向量

    for m = 1:M

        X  = normrnd(mu,ones(n,1));               %n为每个每个X里样本量

        Tx = (mean(X) - 120)/(std(X)/sqrt(n));    %t值

Flag(m) = 1 - (Tx >= -1.96)*(Tx <= 1.96); %在[-1.96 1.96]里面Flag值为0，未落入Flag为1

    end

    Pow(s) = mean(Flag);         %对Flag求均值

end

plot(Mu,Pow);

第5章 

画直方图

1、function Example51(forearm,K)

subplot(1,2,1)

% The hist function optionally returns the

% bin centers and frequencies.

[n,x] = hist(Y,K);

% Plot and use the argument of width=1

% to produce bars that touch.

bar(x,n,1);

axis square

title('Frequency Histogram')

% Now create a relative frequency histogram.

% Divide each box by the total number of points.

subplot(1,2,2)

bar(x,n/140,1)

title('Relative Frequency Histogram')

axis square

2、function Example52(Y,K)

% Get parameter estimates for the normal distribution.

mu = mean(Y);

v = var(Y);

% Obtain normal pdf based on parameter estimates.

xp = linspace(min(Y),max(Y));

yp = normpdf(xp,mu,sqrt(v));

% Get the information needed for a histogram.

[nu,x] = hist(Y,K);

% Get the widths of the bins.

h = x(2)-x(1);

bar(x,nu/(length(Y)*h),'w')

hold on

plot(xp,yp,'r')

hold off

QQ图 Example53

% Generate the random variables.

x = normrnd(0,1,1000,1);

y = trnd(5,1000,1);

% Find the order statistics.

xs = sort(x);

ys = sort(y);

% Now construct the q-q plot.

plot(xs,ys,'o')

xlabel('X - Standard Normal')

ylabel('Y - Standard Normal')

axis equal

hold on 

ezplot('x')

hold off

泊松分布

% Example 5.7

k=0:6;  % vector of counts

n_k = [156 63 29 8 4 1 1];

N=sum(n_k);

% get phi(n_k) for plotting

phik=log(factorial(k).*n_k/N);

% find the counts that are equal to 1

% plot these with the symbol 1

% plot rest with a symbol

plot(k,phik,'o')

% add some whitespace to see better

axis([-0.5 max(k)+1 min(phik)-1 max(phik)+1])

xlabel('Number of Occurrences - k')

ylabel('\\phi (n_k)')

由例题57延伸出来的 

function Example57a(X)

Xmax = max(X);

nk   = zeros(length(Xmax));

for k = 0:Xmax

    nk(k + 1) = sum(X == k);

end

K  = 0:max(X);

Ek = factorial(K).*nk;

id = (Ek ~= 0);

Ek = Ek(id);

phik = log(Ek/sum(nk));

plot(K(id),phik,'o')

% add some whitespace to see better

axis([-0.5 max(K)+1 min(phik)-1 max(phik)+1])

xlabel('Number of Occurrences - k')

ylabel('\\phi (n_k)')

第四章

1、function X= Example42(pp)  %定义函数

pp=[.3;.2;.4;.1];

X0=[100;-100;50;20];

cc = cumsum(pp);  %对pp累加求和

u  = unifrnd(0,1) %定义u为[0，1]上的均匀分布

X=X0(sum((u - cc)>0)+1) %求跳跃点的位置，判断u-cc>0的个数加总为跳跃点,对应取X0中的元素

2、pp图 function [x, xy, irv, rej, rejy, irej] = Example44(c)

c=2;

n = 1000; % Generate 100 random variables.

% Set up the arrays to store variates.

x = zeros(1,n); % random variates

xy = zeros(1,n);% corresponding y values

rej = zeros(1,n);% rejected variates

rejy = zeros(1,n); % corresponding y values

irv = 1;

irej = 1;

while irv <= n

    y = rand(1); % random number from g(y)

    u = rand(1);% random number for comparison

if u <= 2*y/c

        x(irv) = y;

        xy(irv) = u*c;

        irv = irv+1;

    else

        rej(irej) = y;

rejy(irej) = u*c ;% really comparing u*c<=2*y

        irej = irej + 1;

    end

end

plot(x,xy,'o')

hold on

plot(rej,rejy,'x')

hold off

3、function Example44a(X)

X=10;

xx = 0:0.01:1;

F0 = xx.^2;

Fn = zeros(length(xx),1);

for k = 1:length(xx)

    Fn(k) = FunEDF(xx(k),X);

end

plot(F0,Fn,'x')

4、function DemonExample42  %验证Example42函数的正确性

tic;  %计算计算的时间 tic~toc

n = 10000  %n为计算的次数

pp = [.3;.2;.5]

hp = SubDemonExample42(n,pp) 

toc;

function hp = SubDemonExample42(n,pp)

X  = zeros(n,1);

for k = 1:n

    X(k) = SubExample42(pp);

end

X0 = unique(X);  %取X0中不重复的值

hp = zeros(length(X0),1);

for k = 1:length(X0)

    hp(k) = mean(X == X0(k)); %mean为均值，即首先判断X==X0（k），当X=X0（k）时计为1，不等是计为0，之后取10000次[0，1]的均值mean

end

function [X, u] = SubExample42(pp)  %子函数（私有函数，只能在该程序中调用）

cc = cumsum(pp);

u  = unifrnd(0,1);

X = sum((u - cc)>0);

function [X, u] = SubExample42_1(pp) 

pp=[.3;.2;.5];

cc = cumsum(pp);

u  = unifrnd(0,1);

X = sum((u - cc)>0);

第三章 

1、function Fn=Example338(x,X) %定义经验分布函数

Fn=mean(x>=X);

2、%load carsmall

function Example338_1(X)

x=5:.1:50;

Fn=zeros(length(x),1);7

for k=1:length(x)

    Fn(k)=Example338(x(k),X);

end

plot(x,Fn)

3、%load carsmall 汽车行驶的例子，验证函数的分布是否为正态分布

function Example338_2(X)

x=5:.1:50;

hmu=mean(X);

hstd=std(X);

F0=normcdf(x,hmu,hstd);

Fn=zeros(length(x),1);

for k=1:length(x)

    Fn(k)=Example338(x(k),X);

end

plot(x,F0,'r')

plot(x,Fn)

4、function DemoExample338(n) %定义函数名

n=100;

X=normrnd(zeros(n,1),1);

x=-3:.1:3;  %给定x的范围

F0=normcdf(x);

Fn=zeros(length(x),1);    %定义x长度行一列的0矩阵

for k=1:length(x)

    Fn(k)=Example338(x(k),X);

end

plot(x,F0,'r');

hold on 

plot(x,Fn)

第二章

function x = chap_example1(A)

A=[2 -1 .3;-1 .5 .4;3 .4 4];

[E,v] = eig(A);

x     = max(diag(v));

mu=[-2;3];

sig=[1 0.75;0.75 1];

X=mvnrnd(mu,sig,500);

plot(X(:,1)-mu(1),X(:,2)-mu(2),'x');

x=-5:.01:5;

y=0.75*x;

hold on

plot(x,y,'r')

hold off

function hN=example_2_4_1(p0,p)    

p0=0.2

p=0.01

hN=100;

pp=0;

while pp    hN=hN+1;

    pp=1-binocdf(9,hN,p);

end

function [x y]=example1(A)

A=[2 -1 .3;-1 .5 .4;3 .4 4]

[E,v]=eig(A)

x=max(diag(v))

y=min(diag(v))

function [x y]=example2(A)

A=[2 -1 .3;-1 .5 .4;3 .4 4];

v=eig(A)

x=max(v)

y=min(v)

x=0:.1:100;

y1=gamma(x);

plot(x,y1);

y2=gampdf(x,4,2);

plot(x,y2);

id=isnan(MPG); %去空值

id=isnan(MPG)==0;

MPGclean=MPG(id)

%绘制不同均值的正态分布图像，pdf代表概率分布函数（probility distribution function）

mu1=1;

mu2=2;

mu3=3;

x=-10:.1:10;

y1=normpdf(x,mu1,1);

plot(x,y1,'r');

y2=normpdf(x,mu2,1);

y3=normpdf(x,mu3,1);

hold on

plot(x,y2,'g');

hold on

plot(x,y3,'y');

sig1=0.5;

sig2=1;

sig3=2;

mu1=-0.3;

mu2=0;

mu3=0.3;

x=-10:.1:10;

ezplot(@(x)normpdf(x,mu1,sqrt(sig1)));

hold on

ezplot(@(x)normpdf(x,mu2,sqrt(sig2)));

hold on

ezplot(@(x)normpdf(x,mu3,sqrt(sig3)));

function Example44

c=2;

n = 1000; % Generate 100 random variables.

% Set up the arrays to store variates.

x = zeros(1,n); % random variates

xy = zeros(1,n);% corresponding y values

rej = zeros(1,n);% rejected variates

rejy = zeros(1,n); % corresponding y values

irv = 1;

irej = 1;

while irv <= n

    y = rand(1); % random number from g(y)

    u = rand(1);% random number for comparison

if u <= 2*y/c

        x(irv) = y

        xy(irv) = u*c;

        irv = irv+1

    else

        rej(irej) = y;

rejy(irej) = u*c; % really comparing u*c<=2*y

        irej = irej + 1

    end

end

plot(x,xy,'o')

hold on

plot(rej,rejy,'x')

hold off%调用子函数时要删除

[x, xy, irv, rej, rejy, irej] ;

A=x;

X=A;

xx = 0:0.01:1;

F0 = xx.^2;

Fn = zeros(length(xx),1);

for k= 1:length(xx)

    Fn(k) = FunEDF(xx(k),X);

end

figure

plot(F0,Fn,'*')

function Fn = FunEDF(x,X)

Fn = mean(X <= x);