clc;
clear;
close all;
warning off;
addpath(genpath(pwd));
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;
