SigmaPlot Transform to Compute Shelf Life

 

' TRANSFORM TO COMPUTE SHELF LIFE
   'X DATA IS PLACED IN COLUMN X_COL (3 OR MORE DATA POINTS)
   'Y DATA IS PLACED IN COLUMN Y_COL
   'RESULTS ARE PLACED IN COLUMNS RES THROUGH RES+7
      'COLS RES & RES+1 CONTAIN THE REGRESSION LINE
      'COL RES+2 CONTAINS THE LOWER CONFIDENCE LINE
      'COLS RES+3 & RES+4 CONTAIN THE (T90,90) POINT FOR THE DROP LINES
      'COLS RES+5 & RES+6 CONTAIN THE T90 VALUE
      'COL RES+7 IS A WORKING COLUMN
   'COPY A PAGE TEMPLATE ON YOUR DATA TO CREATE THE GRAPH
' INPUT
X_COL=1           'COLUMN NUMBER FOR X DATA
Y_COL=2           'COLUMN NUMBER FOR Y DATA
RES=4              'FIRST RESULTS COLUMN
' PROGRAM
Z=1.96             'Z FOR 95% CONFIDENCE
Y0 = 90            'Y SPECIFICATION LIMIT (%)
X1=COL(X_COL)      'DEFINE X VALUES
Y1=COL(Y_COL)      'DEFINE Y VALUES
   'ROWWISE DELETE MISSING VALUES
FOR ROW = 1 TO SIZE(X1) DO
   CELL(RES+7,ROW)=MISSING(BLOCK(X_COL, ROW, Y_COL, ROW))
END FOR
X=IF(COL(RES+7)=0,X1)
Y=IF(COL(RES+7)=0,Y1)
N=SIZE(X)           'NUMBER OF DATA POINTS
V=N-2                 'N MUST BE > 2
XBAR=MEAN(X)                     'MEAN OF X
DENOM=TOTAL((X-XBAR)^2)          'SUM OF SQS ABOUT MEAN
ALPHA=TOTAL(X^2)/(N*DENOM)       '1,1 COEFF OF (X'X)^-1
BETA=-XBAR/DENOM                 '1,2 COEFF OF (X'X)^-1
DELTA=1/DENOM                    '2,2 COEFF OF (X'X)^-1
R1=TOTAL(Y)                      '1ST ROW OF X'Y
R2=TOTAL(X*Y)                    '2ND ROW OF X'Y
B0=ALPHA*R1+BETA*R2              'INTERCEPT PARAMETER
B1=BETA*R1+DELTA*R2              'SLOPE PARAMETER
   'COMPUTE T VALUE
T123=Z+(Z^3+Z)/(4*V)+(5*Z^5+16*Z^3+3*Z)/(96*V^2)
T4=(3*Z^7+19*Z^5+17*Z^3-15*Z)/(384*V^3)
T5=79*Z^9+776*Z^7+1482*Z^5-1920*Z^3-945*Z
T6=27*Z^11+339*Z^9+930*Z^7-1782*Z^5-765*Z^3+17955*Z
T1=T123+T4+T5/(92160*V^4)+T6/(368640*V^5)
T=IF(V=1, 12.706, IF(V=2, 4.303,T1))
   'ESTIMATE OF S
S=SQRT(TOTAL(((Y-(B0+B1*X))^2))/V)
   'QUADRATIC EQUATION FOR EXPIRATION TIME
DELTA0 = B0-Y0
A = DELTA - (B1/(T*S))^2
B = 2*BETA - 2*B1*DELTA0/(T*S)^2
C = ALPHA - (DELTA0/(T*S))^2
B24AC=B^2-4*A*C
ROOT1=(-B + SQRT(B24AC))/(2*A)
ROOT2=(-B - SQRT(B24AC))/(2*A)
   'FIND APPROPRIATE ROOT
MINX=0
R={ROOT1,ROOT2}
ROOT=IF(S=0, (Y0-B0)/B1, MAX(IF(R<(Y0-B0)/B1,R)))
MAXX=IF(S=0,1.1*ROOT, IF(B1>0 OR B24AC <0, MAX(X), 1.1*ROOT))
XREG=DATA(MINX,MAXX,(MAXX-MINX)/20)
'XREG=DATA(MINX,MAXX,1)  'USED FOR COMPARISON WITH SIGMASTAT
'XREG=X1                 ' USED FOR COMPARISON WITH SAS
YREG=B0+B1*XREG
' CONFIDENCE LIMITS
TERM=ALPHA+2*BETA*XREG+DELTA*XREG^2
CONF_LIM=SQRT(TERM)
LOW_CONF=YREG-T*S*CONF_LIM        ;LOWER LIMIT
' OUTPUT
   'REGRESSION
COL(RES)=XREG           ' X VALUES OF REGRESSION LINE
COL(RES+1)=YREG         ' Y VALUES OF REGRESSION LINE
   'LOWER CONFIDENCE INTERVAL
COL(RES+2)=LOW_CONF     ' LOWER CONFIDENCE LIMIT
'COL(RES+9)=LOW_CONF    ' LOWER CONFIDENCE LIMIT
   'DEGENERATE LINE PLOT FOR SPECIFICATION LIMIT DROP LINES
CELL(RES+3,1)=IF(N<3, "N MUST > 2", IF(B1>0 OR B24AC<0, "NO SOLUTION", ROOT))
CELL(RES+4,1)=Y0
   'SHELF LIFE
CELL(RES+5,1)= "   T90  =  "
CELL(RES+6,1) = IF(B1>0, " + INFINITY",
IF(B24AC<0, " NO SOLUTION", ROOT))

Try SigmaPlot FREE for 30 Days!