本文由 资源共享网 – 08i8 发布，转载请注明出处，如有问题请联系我们！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(t,y,'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(y,2) < m)
   y = y(:, ones(1,m));
end
minyaxis  = log(log(1./(1 - 0.25 ./n)));
maxyaxis  = log(log(1./(0.25 ./n)));
q1x = prctile(t,25);
q3x = prctile(t,75);
q1y = prctile(y,25);
q3y = prctile(y,75);
qx = [q1x; q3x];
qy = [q1y; q3y];
c = zeros(m,2);
mx = [minx maxx];
mx = mx(ones(m,1),:);
my = zeros(m,2);
for k = 1:m
   c(k,:) = polyfit(log(qx(:,k)),qy(:,k),1);
   my(k,:) = polyval(c(k,:),log(mx(k,:)));
end
mx = mx';
my = my';
set(gca,'YTick',tick,'YTickLabel',label,'XScale','log');
set(gca,'YLim',[minyaxis maxyaxis],'XLim',[minxaxis maxxaxis]);
xlabel('Data');
ylabel('Probability');
title('Weibull Probability Plot');
grid on;