diff -rupN scripts/statistics/distributions/wblplot.m scripts-new/statistics/distributions/wblplot.m
--- scripts/statistics/distributions/wblplot.m	1970-01-01 01:00:00.000000000 +0100
+++ scripts-new/statistics/distributions/wblplot.m	2014-11-22 00:36:04.000000000 +0100
@@ -0,0 +1,288 @@
+## 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 <http://www.gnu.org/licenses/>.
+
+## -*- 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 <m95vebj@gmail.com>
+## Created: 2014-11-11
+## 2014-11-22 Updated argin check
+
+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  
+    else
+      error("frequency must be a 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 median rank
+  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 
+
+  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))));
+  poly = polyfit(datx, daty, 1);
+  beta_rry = poly(1);
+  
+  eta_rry = exp(-(poly(2) / beta_rry));
+
+  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;
+  xmaxview=10^ceil(bb);
+  if log10(min(wbdat(:,1)))-floor(log10(min(wbdat(:,1))))<0.5
+    xminview=10^(floor(log10(min(wbdat(:,1))))-1);
+  else
+    xminview=10^floor(log10(min(wbdat(:,1))));      
+  end
+  if xminview<1
+    xminview=1;
+  end
+  if min(wbdat(:,2))>0.20
+    yminview = log(log(1/(1-0.1)));
+  elseif min(wbdat(:,2))>0.03
+    yminview = log(log(1/(1-0.01)));
+  else
+    yminview = log(log(1/(1-0.001)));
+  end
+  xmin=1;
+
+  ymin= log(log(1/(1-0.001)));
+
+
+  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); % log(log(1 / (1 - 10 ^ (-k))));
+      ymajorgrid(k) = line(x,y,'LineStyle','-','Marker','none','Color',[1 0.75 0.75],'LineWidth',0.1);
+    end
+
+    %' Y Minor grids
+    x(1, 1) = xmin;
+    x(2, 1) = xmax*10;
+    for m = 1 : 4
+      for k = 1 : 8
+	%' Y Minor grids - 0.2 - 0.9/Y Minor grids - 0.02 - 0.09/Y Minor grids - 0.002 - 0.009
+	y(1, 1) = log(log(1 / (1 - ((k + 1) / (10 ^ m)))));
+	y(2, 1) = y(1, 1); % log(log(1 / (1 - ((k + 1) / (10 ^ m)))));
+	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 = 0 : log10(xmax)
+      %'X major grids
+      x(1, 1) = 10 ^ m;
+      x(2, 1) =  x(1, 1);% 10 ^ k;
+      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);
+  cb_high=[];
+  cb_low=[];
+  if ~isempty(confint)
+    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];
+    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))));
+    p1=polyfit(log(conf(1:2,1)),confy(1:2,1),1);
+    y1=polyval(p1,log(1));
+    p2=polyfit(log(conf(1:2,1)),confy(1:2,2),1);
+    y2=polyval(p2,log(1));
+    confu=[1 y1 y2;conf(:,1) confy]; # Maybe not totally correct but it looks nicer to extend the lines a bit.
+    h2=plot(confu(:,1),confu(:,2:3),'LineStyle','-','Marker','none','Color',[1 0.25 0.25],'LineWidth',1);
+  end
+  rsq=corr(datx,daty);
+  if showlegend==1
+    s1=sprintf(' RRY\n \\beta=%.3f \n \\eta=%.2f  \n \\rho=%.4f',beta_rry,eta_rry,rsq);
+    if ~isempty(cb_high) &&  ~isempty(cb_low)
+      s2=sprintf('CB_H=%.2f',cb_high);
+      s3=sprintf('CB_L=%.2f',cb_low);
+      legend([h;bestfitrry;h2],"Data",s1,s2,s3,"location","northeastoutside");
+    else
+      legend([h;bestfitrry],"Data",s1,"location","northeastoutside");
+    end
+  end
+  axis([xminview xmaxview yminview log(log(1/(1-0.99)))]);
+  hold off
+  
+  if nargout >= 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
