matlab编的一个用于对统计数据进行威布尔分布的估计，并对其进行三参数的评估与计算77511956Three-parameters-

matlab编的一个用于对统计数据进行威布尔分布的估计，并对其进行三参数的评估与计算。

clc;
clear;
t=[17.8 21.3 23.8 25.9 27.4 29.4 30.6 32.3 33.5 34.9 36.6 38.5 39.7 41.2 43.4 44.5 47 48.8 52.5 61.4];%威布尔书
%%%%%%%%RRY
n=numel（t）;
crep=3000;%迭代次数
r=0;
for j=3:crep
for i=1:n
    f（i）=（i-0.3）/（n+0.4）;
    y（i）=log（log（1/（1-f（i））））;
    x（i）=log（t（i）-r）;
end
for i=1:n
    c1（i）=（x（i）-mean（x））*（y（i）-mean（y））;
    c2（i）=（x（i）-mean（x））^2;
end
c3（j）=sum（c1）/sum（c2）;
c=c3（j）;
d1（j）=exp（-mean（y）/c+mean（x））;
d=d1（j）;
 for i=1:n   
Q_old（i）=（t（i）-r-d*（-log（1-f（i）））^（1/c））.^2;   
r1（i）=（t（i）-d*（-log（1-f（i）））^（1/c））/n;
 end
r2（j）=sum（r1）;
Q（j）=sum（Q_old）;
r=r2（j）;
if Q（j-2）>Q（j-1）&Q（j-1）<=Q（j）  %终止条件
    break
end 
end
canshu=[d c r];%分别为位置参数，形状参数，尺度参数
for i=1:n
Ms_yx（i）=（y（i）-（c*x（i）-c*log（d）））^2;%RRY误差平方
end
Ms_yx1=sum（Ms_yx）;%RRY误差平方和
%%%%%%%%RRX
for j=3:crep
for i=1:n
    f（i）=（i-0.3）/（n+0.4）;
    y（i）=log（log（1/（1-f（i））））;
    x（i）=log（t（i）-r）;
end
for i=1:n
    c1（i）=（y（i）-mean（y））^2;
    c2（i）=（x（i）-mean（x））*（y（i）-mean（y））;
end
c3（j）=sum（c1）/sum（c2）;
c=c3（j）;
d1（j）=exp（mean（x）-mean（y）/c）;
d=d1（j）;
 for i=1:n   
Q_old（i）=（t（i）-r-d*（-log（1-f（i）））^（1/c））.^2;   
r1（i）=（t（i）-d*（-log（1-f（i）））^（1/c））/n;
 end
r2（j）=sum（r1）;
Q（j）=sum（Q_old）;
r=r2（j）;
if Q（j-2）>Q（j-1）&Q（j-1）<=Q（j）  %终止条件
    break
end
end
canshu1=[d c r];
for i=1:n
Ms_xy（i）=（x（i）-（log（d）+y（i）/c））^2;%RRX误差平方
end
Ms_xy1=sum（Ms_xy）;%RRX误差平方和
%%%%%%%两种方法的误差平方和比较，得到最优参数
if Ms_xy1>Ms_yx1
    true=canshu;
else
    true=canshu1;   
end
fprintf（‘\n位置参数为：%f‘true（1））;
fprintf（‘\n形状参数为：%f‘true（2））;
fprintf（‘\n尺度参数为：%f‘true（3））;
%%%%%%%%%画图
figure;
plot（ty‘b*‘）
hold on
plot（t（-true（2）*log（true（1））+true（2）.*log（t-true（3）））‘k-‘）
%关于减去位置参数后与Y的关系图
plot（（t-true（3））y‘c*‘）
hold on
plot（（t-true（3））（-true（2）*log（true（1））+true（2）.*log（t-true（3）））‘m-‘）
legend（‘失效数据‘‘拟合曲线‘‘失效数据-位置参数‘‘（失效数据-位置参数）拟合的曲线‘）
%%%%%%weibull概率纸
[n m] = size（t）;
if n == 1
   t = t‘;
   n = m;
   m = 1;
end
[sx i]= sort（t）;
minx  = min（sx（:））;
maxx  = max（sx（:））;
range = maxx-minx;
if range > 0
  minxaxis  = 0;
  maxxaxis  = maxx+0.025*range;
else
  minxaxis  = minx - 1;
  maxxaxis  = maxx + 1;
end
p     = [0.001 0.003 0.01 0.02 0.05 0.10 0.25 0.5...
         0.75 0.90 0.96 0.99 0.999];
label1= str2mat（‘0.001‘‘0.003‘ ‘0.01‘‘0.02‘‘0.05‘‘0.10‘‘0.25‘‘0.50‘）;
label2= str2mat（‘0.75‘‘0.90‘‘0.96‘‘0.99‘ ‘0.999‘）;
label = [label1;label2];
tick  = log（log（1./（1-p）））;
%%%%% Y坐标的设定
if （~any（isnan（t（:））））
   eprob = [0.5./n:1./n:（n - 0.5）./n]‘;
else
   nvec = sum（~isnan（t））;
   eprob = repmat（（1:n）‘ 1 m）;
   eprob = （eprob-.5） ./ repmat（nvec n 1）;
   eprob（isnan（sx）） = NaN;
   n = max（nvec）;  % sample size for setting axis limits
end
y  = log（log（1./（1-eprob）））;
if （size（y2） < m）
   y = y（: ones（1m））;
end
minyaxis  = log（log（1./（1 - 0.25 ./n）））;
maxyaxis  = log（log（1./（0.25 ./n）））;
q1x = prctile（t25）;
q3x = prctile（t75）;
q1y = prctile（y25）;
q3y = prctile（y75）;
qx = [q1x; q3x];
qy

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3574  2011-05-16 15:47  Three parameters complete data about the least square method.m

----------- ---------  ---------- -----  ----

                 3574                    1