Lognormal Plotting of Single Data (Example Problem 5.1), MATLAB

%% Lognormal Plotting (Single) Example
% Coded by Reuel Smith 2015-2017
% v. MATLAB R2015b through 2017a
% ========================================================================
% Example Problem 5.1
% Five specimens of a new corrosion-resistant material are tested for 240
% hours in a highly corrosive environment.  The density of the material is
% 7.6 g/cm3, and the exposed surface area of each specimen is 4.3 cm2.  At
% the end of the test period, the measured weight losses in mg were 11.1,
% 10.4, 12.1, 11.4, and 9.8 (assume no measurement error).  If a
% degradation of 1 mm or more results in a structural failure, predict the
% failure times for the five specimens.  Based on these failure times,
% determine the probability distribution that best represents the life of
% the material.
%
% Solution:
% TTF = 70659, 75415, 64820, 68800, and 80033
% The proposed probability solutions for this problem are Weibull and
% Lognormal.  This routine will process the Lognormal solution for the
% parameters log mean and log sigma, and the Lognormal fit probability plot
% ========================================================================
clc
clear
% ========================================================================
% PART 1: Initialize data and constants
% ========================================================================
% sorting the times to failure in ascending order
TTF = sort([70659 75415 64820 68800 80033])';

% set the length of the TTF vector or the number of failures
n = length(TTF);

% CDF or Failure Percent for probability plot under the Kimball probability
% plotting position: (i - 0.375)/(n + 0.25).  Other plotting positions may
% be assumed:
% Midpoint - (i - 0.5)/n
% Mean - i/(n + 1)
% Median - (i - 0.3)/(n + 0.4)
i = 1:n;                            % uncensored data points from 1 to n
F = 100.*((i'-0.375)./(n+0.25));
% Lognormal scale conversion of the Failure Percent data
FB = icdf('Normal',F./100,0,1);

% reliability for probability plot
R = (100-F)./100;

% Scale conversion to Lognormal plot scale
% Upper and lower bounds of the Percent Failure axis in percent
fc = [.1 99.9];
% Lognormal scale conversion of the upper and lower bounds of the Percent 
% Failure axis 
fcB = icdf('Normal',fc./100,0,1);
% Percent Failure ticks cited between the lower and upper bounds
Pticks1 = [0.1:0.1:1 2:1:10 20:10:90 95 99 99.9];
% Lognormal scale conversion of Percent Failure ticks cited between the
% lower and upper bounds 
Pticks = icdf('Normal',Pticks1./100,0,1);
% ========================================================================
% PART 2: Solve for Lognormal Parameters log mean and log sigma
% ========================================================================
% Use a polynomial fit (polyfit) to draw the best fit line (fc) through the
% known times to failure TTF from F(t) = 1% to F(t) = 99%
% polynomial fit for LN(TTF) and F
pb2 = polyfit(log(TTF),F,1);
% the log mean of the data
meant = (50-pb2(2))/pb2(1);
% the 84% value for LN(t)
t84 = (84-pb2(2))/pb2(1);
% the log-sigma
sigmat = t84-meant;
lognparam = lognfit(TTF);
disp(['Lognormal parameters: log mean = ',num2str(meant),' and log sigma = ',num2str(sigmat)])
% ========================================================================
% PART 3: Lognormal probability plot
% ========================================================================
% Plot points at 0.1% and 99.9%
% These have a different polynomial pair than that which was used to find
% the meant and sigmat
pb2_B = polyfit([meant meant+sigmat],[icdf('Normal',0.5,0,1) icdf('Normal',0.84,0,1)],1);
ttfc = exp([(icdf('Normal',0.001,0,1)-pb2_B(2))/pb2_B(1) (icdf('Normal',0.999,0,1)-pb2_B(2))/pb2_B(1)]);
signs1 = [floor(log10(ttfc(1))):1:ceil(log10(ttfc(2)))];

figure(1)
semilogx(TTF,FB,'r*',ttfc,fcB,'r-');
xlabel('Time (hours)')
ylabel('Percent Failure')
grid on
set(gca,'YTick',Pticks,'YTickLabel',{'0.10';'0.20';' ';' ';'0.50';' ';' ';' ';' ';'1.00';'2.00';' ';' ';'5.00';' ';' ';' ';' ';'10.00';'20.00';'30.00';'40.00';'50.00';'60.00';'70.00';'80.00';'90.00';'95.00';'99.00';'99.90'})
set(gca,'XTickLabel',{num2str(10.^signs1')})
legend('Data','Lognormal best fit line','Location','Northwest')
xlim([10^signs1(1) 10^signs1(end)])
ylim([icdf('Normal',0.001,0,1) icdf('Normal',0.999,0,1)])
set(gcf,'color','w');