中国DOS联盟

%打印T分布表 {19:32 2008-5-13|plp626}
function m() %定义主函数
!set/p=%time%<nul>time
v=[0.25 0.10 0.05 0.025 0.01 0.005];%设置分位数向量
fprintf('n╲α')
for i=1:1:6
fprintf('%10.4f',v(i))
end

fprintf('\n')
for j=1:1:45                       %设置自由度范围1-45
fprintf('%5.1d  ',j)
   for i=1:1:6
      fprintf('%10.4f',dt(v(i),j))
      if i==6 fprintf('\n')
      end
   end
end
!echo\ %time%>>time
!for /f "tokens=*" %a in (time)do timediff %a 0

%对自由度为n,求积分值为v（分为数）时的积分下限值。
function r=dt(v,n)   %%//////dt(v,n)
h=1;                 %初始步长
a=0.5;               %初值a，然后迭代到我们要求的值。
ta=t(a,n);           %初值的积分值:0.5为积分下限。
%%%%这里采用变步长3-4-5逼近法（原创）
while abs(ta-v)>0.000001 %设置精度
        while ta>v
                a=a+h;ta=t(a,n);
        end
        h=0.3*h;a=a-h;
        while ta<v
                a=a-h;ta=t(a,n);
        end
        h=0.4*h;a=a+0.5*h;
end
r=a;

%计算自由度为n时[a0,inf]的t积分。对于指定的n,t(a,n)为递减函数。
function r=t(a0,n)                %%/////t(a0,n)
syms x
f=gamma((n+1)/2)/(sqrt(n*pi)*gamma(n/2))*(1+x*x/n)^(-(n+1)/2);%T分布密度函数

f=simplify(f);
g=int(f,x,a0,inf);
r=numeric(g);

%打印正态分布表  
%$作者：****** 时间：2008-5-23 $
function f()
p=0:9;                                        %第2位小数
fprintf('   z  ')
for i=1:10
fprintf('%-8.1d',p(i))
end
fprintf('\n')
for z=0:0.1:2.0                               %第一位小数
   fprintf('%4.1f',z)
   for i=1:1:10
      fprintf('%8.4f',int_f(z+p(i)/100))
      if i==10 fprintf('\n')
      end
   end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r=int_f(a)
syms x
r=double(int(1/sqrt(2*pi)*exp(x^2/-2),-inf,a));