1、第三次上机作业计量上机实验报告（3）计量经济学上机实验报告实验内容：任务1-5实验工具：Matlab2008实验目的：学会使用if语句、循环、函数、积分和极大似然估计的基本方法实验过程：任务1:求1到10000中能被3整除的数的和.Matlab代码s=0 for i=1:1:10000 p(i)=i; q(i)=mod(i,3); if q(i)=fix(q(i) s=s+p(i); end; end;s结果s = 16668333任务2：用function()命令编写正态分布的密度函数，并分别计算N(0,1)，N(0,2)，N(0,3)小于1部分的面积。Matlab 代码 y2=normcd
2、f(x,0,1); plot(x,y1,r-,x,y2,g-); plot(x,y1,r-,x,y2,g-);legend(密度函数,分布函数);J=int(normpdf(x,0,1),x,-inf,1);syms x;J1=int(normpdf(x,0,2),x,-inf,1);VJ1=vpa(J1); syms x;J2=int(normpdf(x,0,3),x,-inf,1);VJ2=vpa(J2); VJ=vpa(J) VJ1 VJ2结果VJ = 0.84134474606854303612337332991661VJ1 =0.691462461274013175581267062
3、14263VJ2 =0.63055865981823639009615938887623概率密度函数曲线图和概率分布曲线图任务3：生成10000个标准正态分布的随机数，分别计算大于1.96，小于-1.96和-1.96到1.96之间的数各有多少个。Matlab 代码s2=0;for i=1:1:10000 if(x3(i)1.96) s2=s2+1; end;end; s2结果s2 = 245任务4：生成y=500+125*x+epslon，epslonN(0,10)的1000个数据点，并用极大似然估计进行回归。Matlab代码Command 主函数epslon=normrnd(0,10,100
4、0,1);x=randn(1000,1);y=500+125*x+epslon;para,standard_deviation=my_mle(mynormpdfsum,500;125;10,y,x);z=para(1)+para(2)*x;plot(x,y,o);hold on;plot(x,z,-);Define mynormpdfsumfunction f=mynormpdfsum(para,num,y,x)yy=1/sqrt(2*pi)/para(3)*exp(-(y-para(1)-para(2)*x).2/2/para(3)2);if num=1 %(note: it must be
5、 set to 1) f=log(yy);else f=-sum(log(yy);endDefine my_mlefunction para,standard_deviation,fv=my_mle(fun,para0,varargin)%estimate parameters and standard errors when using maximium likelihood%estimation(MLE)%input%fun: a function defined by users for calculating log probability density%function (pdf)
6、 and negative sum of logarithm of pdf%para0 : given initial parameters% varargin: other needed parameters required by fun%output%para: estimated parameters% standard_deviation: standard deviations of estimated parameters% fv: maximized likelihood function value%para0=para0(:);para,fv=fminsearch(fun,
7、para0,2,varargin:);fv=-fv;n=length(para);d=numericalfirstderivative(fun,para,1,varargin:);standard_deviation=sqrt(diag(pinv(d*d);%function f=numericalfirstderivative(fun,parameter,varargin)% input:% fun: the name of a function% parameter: given parameter with respect to which first-order derivative%
8、is calculated% varargin: other needed parameters required by fun% output:% f: numerical first order derivative of fun at parametern=length(parameter);for i=1:na=zeros(n,1);a(i)=min(parameter(i)*1e-6,1e-5);y1(:,i)=feval(fun,parameter+a,varargin:);y2(:,i)=feval(fun,parameter-a,varargin:);f(:,i)=(y1(:,
9、i)-y2(:,i)/2/a(i);end结果para = 499.4534 125.1685 10.1249standard_deviation = 0.3204 0.3244 0.2104Figure分析上述估计表明：a、b的参数估计分别为499.4534,125.1685,epslon的均值和方差的估计值分别为0.2104、10.1249，这与给出的标准值之间有一定的差值。任务5：生成y=2+ex+epslon，epslonN(0,1)的1000个数据点，并用极大似然估计进行回归。Matlab代码Command 主函数epslon=normrnd(0,10,1000,1);x=randn
10、(1000,1);y=2+exp(x)+epslon;para,standard_deviation=my_mle(mynormpdfsum,2;exp(1);1,y,x);z=para(1)+para(2).x;plot(x,y,o);hold on;plot(x,z,-);title(极大似然估计)Define mynormpdfsumfunction f=mynormpdfsum(para,num,y,x)yy=1/sqrt(2*pi)/para(3)*exp(-(y-para(1)-para(2).x).2/2/para(3)2);if num=1 %(note: it must be
11、 set to 1) f=log(yy);else f=-sum(log(yy);endDefine my_mlefunction para,standard_deviation,fv=my_mle(fun,para0,varargin)%estimate parameters and standard errors when using maximium likelihood%estimation(MLE)%input%fun: a function defined by users for calculating log probability density%function (pdf)
12、 and negative sum of logarithm of pdf%para0 : given initial parameters% varargin: other needed parameters required by fun%output%para: estimated parameters% standard_deviation: standard deviations of estimated parameters% fv: maximized likelihood function value%para0=para0(:);para,fv=fminsearch(fun,
13、para0,2,varargin:);fv=-fv;n=length(para);d=numericalfirstderivative(fun,para,1,varargin:);standard_deviation=sqrt(diag(pinv(d*d);%function f=numericalfirstderivative(fun,parameter,varargin)% input:% fun: the name of a function% parameter: given parameter with respect to which first-order derivative%
14、is calculated% varargin: other needed parameters required by fun% output:% f: numerical first order derivative of fun at parametern=length(parameter);for i=1:na=zeros(n,1);a(i)=min(parameter(i)*1e-6,1e-5);y1(:,i)=feval(fun,parameter+a,varargin:);y2(:,i)=feval(fun,parameter-a,varargin:);f(:,i)=(y1(:,i)-y2(:,i)/2/a(i);endFigure结果para = 1.9873 2.7273 1.0242standard_deviation = 0.0340 0.0138 0.0230分析上述估计表明：a、e的参数估计分别为1.9873，2.7273,epslon的均值和方差的估计值分别为0.0230、1.0242，这与给出的标准值之间有一定的差值。
