第三次上机作业.docx
- 文档编号:30208363
- 上传时间:2023-08-07
- 格式:DOCX
- 页数:9
- 大小:106.29KB
第三次上机作业.docx
《第三次上机作业.docx》由会员分享,可在线阅读,更多相关《第三次上机作业.docx(9页珍藏版)》请在冰豆网上搜索。
第三次上机作业
计
量
上
机
实
验
报
告
(3)
计量经济学上机实验报告
实验内容:
任务1-5
实验工具:
Matlab2008
实验目的:
学会使用if语句、循环、函数、积分和极大似然估计的基本方法
实验过程:
任务1:
求1到10000中能被3整除的数的和.
Matlab代码
s=0
fori=1:
1:
10000
p(i)=i;
q(i)=mod(i,3);
ifq(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=normcdf(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);
symsx;
J1=int(normpdf(x,0,2),x,-inf,1);
VJ1=vpa(J1);
>>symsx;
J2=int(normpdf(x,0,3),x,-inf,1);
VJ2=vpa(J2);
>>VJ=vpa(J)
>>VJ1
>>VJ2
结果
VJ=
0.84134474606854303612337332991661
VJ1=
0.69146246127401317558126706214263
VJ2=
0.63055865981823639009615938887623
概率密度函数曲线图和概率分布曲线图
任务3:
生成10000个标准正态分布的随机数,分别计算大于1.96,小于-1.96和-1.96到1.96之间的数各有多少个。
Matlab代码
s2=0;
fori=1:
1:
10000
if(x3(i)>1.96)
s2=s2+1;
end;
end;
>>s2
结果
s2=
245
任务4:
生成y=500+125*x+epslon,epslon~N(0,10)的1000个数据点,并用极大似然估计进行回归。
Matlab代码
Command主函数
epslon=normrnd(0,10,1000,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');
holdon;
plot(x,z,'--');
Definemynormpdfsum
functionf=mynormpdfsum(para,num,y,x)
yy=1/sqrt(2*pi)/para(3)*exp(-(y-para
(1)-para
(2)*x).^2/2/para(3)^2);
ifnum==1%(note:
itmustbesetto1)
f=log(yy);
else
f=-sum(log(yy));
end
Definemy_mle
function[para,standard_deviation,fv]=my_mle(fun,para0,varargin)
%estimateparametersandstandarderrorswhenusingmaximiumlikelihood
%estimation(MLE)
%input
%fun:
afunctiondefinedbyusersforcalculatinglogprobabilitydensity
%function(pdf)andnegativesumoflogarithmofpdf
%para0:
giveninitialparameters
%varargin:
otherneededparametersrequiredbyfun
%output
%para:
estimatedparameters
%standard_deviation:
standarddeviationsofestimatedparameters
%fv:
maximizedlikelihoodfunctionvalue
%%%%%%%%%%%
%
para0=para0(:
);
[para,fv]=fminsearch(fun,para0,[],2,varargin{:
});
fv=-fv;
n=length(para);
d=numericalfirstderivative(fun,para,1,varargin{:
});
standard_deviation=sqrt(diag(pinv(d'*d)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionf=numericalfirstderivative(fun,parameter,varargin)
%input:
%fun:
thenameofafunction
%parameter:
givenparameterwithrespecttowhichfirst-orderderivative
%iscalculated
%varargin:
otherneededparametersrequiredbyfun
%output:
%f:
numericalfirstorderderivativeoffunatparameter
n=length(parameter);
fori=1:
n
a=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);
end
结果
para=
499.4534
125.1685
10.1249
standard_deviation=
0.3204
0.3244
0.2104
Figure
分析
上述估计表明:
a、b的参数估计分别为499.4534,125.1685,epslon的均值和方差的估计值分别为0.2104、10.1249,这与给出的标准值之间有一定的差值。
任务5:
生成y=2+e^x+epslon,epslon~N(0,1)的1000个数据点,并用极大似然估计进行回归。
Matlab代码
Command主函数
epslon=normrnd(0,10,1000,1);
x=randn(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');
holdon;
plot(x,z,'--');
title('极大似然估计')
Definemynormpdfsum
functionf=mynormpdfsum(para,num,y,x)
yy=1/sqrt(2*pi)/para(3)*exp(-(y-para
(1)-para
(2).^x).^2/2/para(3)^2);
ifnum==1%(note:
itmustbesetto1)
f=log(yy);
else
f=-sum(log(yy));
end
Definemy_mle
function[para,standard_deviation,fv]=my_mle(fun,para0,varargin)
%estimateparametersandstandarderrorswhenusingmaximiumlikelihood
%estimation(MLE)
%input
%fun:
afunctiondefinedbyusersforcalculatinglogprobabilitydensity
%function(pdf)andnegativesumoflogarithmofpdf
%para0:
giveninitialparameters
%varargin:
otherneededparametersrequiredbyfun
%output
%para:
estimatedparameters
%standard_deviation:
standarddeviationsofestimatedparameters
%fv:
maximizedlikelihoodfunctionvalue
%%%%%%%%%%%
%
para0=para0(:
);
[para,fv]=fminsearch(fun,para0,[],2,varargin{:
});
fv=-fv;
n=length(para);
d=numericalfirstderivative(fun,para,1,varargin{:
});
standard_deviation=sqrt(diag(pinv(d'*d)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
functionf=numericalfirstderivative(fun,parameter,varargin)
%input:
%fun:
thenameofafunction
%parameter:
givenparameterwithrespecttowhichfirst-orderderivative
%iscalculated
%varargin:
otherneededparametersrequiredbyfun
%output:
%f:
numericalfirstorderderivativeoffunatparameter
n=length(parameter);
fori=1:
n
a=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);
end
Figure
结果
para=
1.9873
2.7273
1.0242
standard_deviation=
0.0340
0.0138
0.0230
分析
上述估计表明:
a、e的参数估计分别为1.9873,2.7273,epslon的均值和方差的估计值分别为0.0230、1.0242,这与给出的标准值之间有一定的差值。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第三次 上机 作业