## Copyright (C) 2014 Björn Vennberg
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program. If not, see .
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{handle} @var{result}] =} wblplot (@var{data} , @var{censor},
## @var{frequ}, @var{confint}, @var{fancygrid}, @var{showtext})
##
## @deftypefnx {Function File} {[@var{handle} @var{result}] =} wblplot (@var{data},...)
##
##
##
## @noindent
## Plot singel column vector @var{data} on a Weibull probability plot using Rank Regression on Y.
##
## @var{censor} optional parameter is a singel vector of same size as @var{data}
## with 1 for right censored data and 0 for exact observation.
##
## @var{frequ} optional vector same size as data with the number of occurencies for corresponding data.
##
## @var{confint} optional confidens limits ploting upper and lower
## confidens band using beta binomial confidence bounds. If a single
## value is given this will be used such as LOW = a and HIGH = 1 - a.
##
## @var{fancygrid} optional paramter which if set to anything but 1 will turn of the the fancy gridlines.
##
## @var{showlegend} optional paramter that when set to zero(0) turns off the legend.
##
## If one output argument is given a @var{handle} for the data marker and plotlines are returned
## which can be used for further modification of line-, marker-style.
##
## If a second output argument is specified a @var{result} vector with scale,
## shape and correlation factor is returned.
##
## @seealso{}
## @end deftypefn
## Author: Björn Vennberg
## Created: 2014-11-11
## 2014-11-22 Updated argin check
## 2014-11-22 Code clean up, error checking
function [handle result] = wblplot (data , censor=[], frequ=[], confint=[],fancygrid=1,showlegend=1)
[mm nn] = size(data);
if mm > 1 && nn > 1
error ("wblplot can only handle a single data vector")
elseif mm == 1 && nn > 1
data=data(:);
mm = nn;
end
if isempty(frequ)
frequ = ones(mm,1);
N = mm;
else
[mmf nnf]=size(frequ);
if (mmf == mm && nnf == 1) || (mmf == 1 && nnf == mm)
frequ = frequ(:);
N=sum(frequ); % Total number of samples
if any(frequ<=0)
error("frequency vector must be all positive non zero integers")
end
else
error("frequency must be vector of same length as data")
end
end
if isempty(censor)
censor = zeros(mm,1);
else
[mmc nnc]=size(censor);
if (mmc == mm && nnc == 1) || (mmc == 1 && nnc == mm)
censor = censor(:);
else
error("censor must be a vector of same length as data")
end
end
## Determin the order number
wbdat=zeros(length(find(censor==0)),3);
Op = 0;
Oi = 0;
c = N;
nf = 0;
for k = 1 : mm
if censor(k, 1) == 0
nf = nf + 1;
wbdat(nf, 1) = data(k, 1);
for s = 1 : frequ(k, 1);
Oi = Op + ((N + 1) - Op) / (1 + c);
Op = Oi;
c = c - 1;
end
wbdat(nf, 3) = Oi;
else
c = c - frequ(k, 1);
endif
end
# Compute median rank
a=wbdat(:, 3)./(N-wbdat(:, 3)+1);
f=finv(0.5,2*(N-wbdat(:, 3)+1),2*wbdat(:, 3));
wbdat(:, 2) = a./(f+a);
datx = log(wbdat(:,1));
daty = log(log(1 ./ (1 - wbdat(:,2))));
# Rank regression
poly = polyfit(datx, daty, 1);
# Shape factor
beta_rry = poly(1);
# Scale factor
eta_rry = exp(-(poly(2) / beta_rry));
# Determin min-max values of view port
aa=ceil(log10(max(wbdat(:,1))));
bb=log10(max(wbdat(:,1)));
if aa-bb<0.2
aa=ceil(log10(max(wbdat(:,1))))+1;
end
xmax= 10^aa;
if log10(min(wbdat(:,1)))-floor(log10(min(wbdat(:,1))))<0.2
xmin=10^(floor(log10(min(wbdat(:,1))))-1);
else
xmin=10^floor(log10(min(wbdat(:,1))));
end
if min(wbdat(:,2))>0.20
ymin = log(log(1/(1-0.1)));
elseif min(wbdat(:,2))>0.02
ymin = log(log(1/(1-0.01)));
elseif min(wbdat(:,2))>0.002
ymin = log(log(1/(1-0.001)));
else
ymin = log(log(1/(1-0.0001)));
end
ymax= log(log(1/(1-0.999)));
x=[0;0];
y=[0;0];
label = char('0.10','1.00','10.00','99.00');
prob = [0.001 0.01 0.1 0.99];
tick = log(log(1./(1-prob)));
xbf = [xmin;xmax];
ybf = polyval(poly, log(xbf));
newplot();
x(1, 1) = xmin;
x(2, 1) = xmax;
if fancygrid==1
for k = 1 : 4
%' Y major grids
x(1, 1) = xmin;
x(2, 1) = xmax*10;
y(1, 1) = log(log(1 / (1 - 10 ^ (-k))));
y(2, 1) = y(1, 1);
ymajorgrid(k) = line(x,y,'LineStyle','-','Marker','none','Color',[1 0.75 0.75],'LineWidth',0.1);
end
%' Y Minor grids 2 - 9
x(1, 1) = xmin;
x(2, 1) = xmax*10;
for m = 1 : 4
for k = 1 : 8
y(1, 1) = log(log(1 / (1 - ((k + 1) / (10 ^ m)))));
y(2, 1) = y(1, 1);
yminorgrid(k) = line(x,y,'LineStyle','-','Marker','none','Color',[0.75 1 0.75],'LineWidth',0.1);
end
end
#'X-axis grid
y(1, 1) = ymin;
y(2, 1) = ymax;
for m = log10(xmin) : log10(xmax)
x(1, 1) = 10 ^ m;
x(2, 1) = x(1, 1);
y(1, 1) = ymin; %
y(2, 1) = ymax; %
xmajorgrid(k) = line(x,y,'LineStyle','-','Marker','none','Color',[1 0.75 0.75]);
for k = 1 : 8
%' X Minor grids - 2 - 9
x(1, 1) = (k + 1) * (10 ^ m);
x(2, 1) = (k + 1) * (10 ^ m);
xminorgrid(k) = line(x,y,'LineStyle','-','Marker','none','Color',[0.75 1 0.75],'LineWidth',0.1);
end
end
end
set(gca,'XScale','log');
set(gca,'YTick',tick,'YTickLabel',label);
xlabel('Data','FontSize',12);
ylabel('Unreliability, F(t)=1-R(t)','FontSize',12);
title('Weibull Probability Plot','FontSize',12);
set(gcf,'Color',[0.9,0.9,0.9])
set(gcf,'name','WblPlot')
hold on
h=plot(wbdat(:,1),daty,'o');
set(h,'markerfacecolor',[0,0,1])
set(h,'markersize',8)
bestfitrry = line(xbf,ybf,'LineStyle','-','Marker','none','Color',[0.25 0.25 1],'LineWidth',1);
# If requested plot beta binomial confidens bounds
if ~isempty(confint)
cb_high=[];
cb_low=[];
if length(confint)==1
if confint >0.5
cb_high=confint;
cb_low=1-confint;
else
cb_high=1-confint;
cb_low=confint;
end
else
cb_high=confint(2);
cb_low=confint(1);
end
conf=zeros(N+4,3);
betainv = 1 / beta_rry;
N2 = [1:N]';
N2=[0.3;0.7;N2;N2(end)+0.5;N2(end)+0.8]; # Extend the ends a bit
ypos = medianranks(0.5, N, N2);
conf(:, 1) = eta_rry * log(1./ (1 - ypos)).^ betainv;
conf(:, 2) = medianranks(cb_low, N, N2);
conf(:, 3) = medianranks(cb_high, N, N2);
confy=log(log(1./(1-conf(:,2:3))));
confu=[conf(:,1) confy];
if conf(1,1)>xmin # Not totally correct but it looks better to extend the lines a bit.
p1=polyfit(log(conf(1:2,1)),confy(1:2,1),1);
y1=polyval(p1,log(xmin));
p2=polyfit(log(conf(1:2,1)),confy(1:2,2),1);
y2=polyval(p2,log(xmin));
confu=[xmin y1 y2;confu];
end
if conf(end,1)= 2
result = [eta_rry beta_rry rsq];
if ~isempty(confint)
handle = [h; bestfitrry; h2];
else
handle = [h; bestfitrry];
end
end
if nargout == 1
if ~isempty(confint)
handle = [h; bestfitrry; h2];
else
handle = [h; bestfitrry];
end
end
endfunction
function [ ret ] = medianranks (alpha, n, ii)
a=ii./(n-ii.+1);
f=finv(alpha,2*(n-ii.+1),2*ii);
ret=a./(f.+a);
endfunction