(* Content-type: application/vnd.wolfram.mathematica *)

(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)

(* CreatedBy='Mathematica 9.0' *)

(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[       157,          7]
NotebookDataLength[     20274,       1131]
NotebookOptionsPosition[     20026,       1118]
NotebookOutlinePosition[     20370,       1133]
CellTagsIndexPosition[     20327,       1130]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{
Cell["\<\
(* PROBIT2-09 - This program is used to analyze serial time-mortality data 

from bioassay studies. The program reads input from the file PR2input. The 

user should make no changes to the program below. *)



(* This program was developed by:



   Jim Throne

   USDA, ARS, SJVASC

   9611 South Riverbend Avenue

   Parlier CA 93648

   (559) 596-2702

   james.throne@ars.usda.gov



   Please notify me of problems or suggestions.  3 June 2013 *)



Off[]



(* Read the input data set. *)



data = ReadList[\"PR2input\"];



filename = data[[1]];

title1 = data[[2]];

title2 = data[[3]];

notimes = data[[4]];

contdead = data[[5]];

expdead = data[[6]];

intimes = data[[7]];

ncont = data[[8]];

nexp = data[[9]];

graphs = data[[10]];



(* The next line opens the output file for writing. *)



out=OpenWrite[filename,FormatType->OutputForm,PageWidth->78];



(* The next lines put descriptions on the output file. *)



WriteString[out,title1]

WriteString[out,\" \\n\"]

WriteString[out,title2]

WriteString[out,\" \\n\"]

WriteString[out,\" \\n\\n\"]



(* Descriptions of input data.



   notimes   = number of time intervals for which you will enter data

   contdead  = the observed cumulative number of insects dead in the

               control at each time interval, enclosed in brackets

   expdead   = the observed cumulative number of insects dead in the

               treatment at each time interval, enclosed in brackets

   times     = the times at which observations were made, enclosed in

               brackets

   ncont     = the initial number of insects in the control

   nexp      = the initial number of insects in the treatment



  graphs     = decision to print out graphs showing observed and

               predicted probits - yes or no in quotation marks *)





(* The user should not make any changes to the program below. *)



WriteString[out,\"Input data:\\n\\n\"];

WriteString[out,\"Number of time intervals = \"];

Write[out,notimes];

WriteString[out,\"Number dead in controls = \"];

Write[out,contdead];

WriteString[out,\"Number dead in treatments = \"];

Write[out,expdead];

WriteString[out,\"Time intervals = \"];

Write[out,intimes];

WriteString[out,\"Number of controls = \"];

Write[out,ncont];

WriteString[out,\"Number treated = \"];

Write[out,nexp];

WriteString[out,\"\\n\"];



ndist=NormalDistribution[0,1];



loop=1;

Do[Label[loopstart];

If[loop==1,tranx=\"no\",Goto[loop2]];

trany=\"probit\";

WriteString[out,\"\\n\\n*** Probit Transformation ***\\n\\n\\n\"];

Goto[loopend];

Label[loop2];

If[loop==2,tranx=\"no\",Goto[loop3]];

trany=\"logit\";

WriteString[out,\"\\n\\n*** Logit Transformation ***\\n\\n\\n\"];

Goto[loopend];

Label[loop3];

If[loop==3,tranx=\"no\",Goto[loop4]];

trany=\"cll\";

WriteString[out,\"\\n\\n*** CLL Transformation ***\\n\\n\\n\"];

Goto[loopend];

Label[loop4];

If[loop==4,tranx=\"yes\",Goto[loop5]];

trany=\"probit\";

WriteString[out,\"\\n\\n*** Log-Probit Transformation ***\\n\\n\\n\"];

Goto[loopend];

Label[loop5];

If[loop==5,tranx=\"yes\",Goto[loop6]];

trany=\"logit\";

WriteString[out,\"\\n\\n*** Log-Logit Transformation ***\\n\\n\\n\"];

WriteString[out,\" \\n\"];

Goto[loopend];

Label[loop6];

tranx=\"yes\";

trany=\"cll\";

WriteString[out,\"\\n\\n*** Log-CLL Transformation ***\\n\\n\\n\"];

WriteString[out,\" \\n\"];

Label[loopend];



(* Delete data for time periods when no mortality occurred or in

which all insects were dead.  Apply log10 transformation to time. *)



If[tranx==\"yes\",times=Log[10,intimes] //N,times=intimes //N];

input=Table[0,{i,notimes},{j,3}];

input[[1,1]]=times[[1]];

input[[1,2]]=contdead[[1]];

input[[1,3]]=expdead[[1]];

inc=1;

Do[If[expdead[[i]]==expdead[[i-1]],Goto[elim],Goto[add]];

Label[add];

If[expdead[[i]]==nexp,Goto[elim],inc=inc+1];

input[[inc,1]]=times[[i]];

input[[inc,2]]=contdead[[i]];

input[[inc,3]]=expdead[[i]];

Goto[end1];

Label[elim];

If[expdead[[i]]==0,input[[1,1]]=times[[i]]];

If[expdead[[i]]==0,input[[1,3]]=contdead[[i]]];

Label[end1];

,{i,2,notimes}];



(* Set up vectors containing the time intervals (t), the number dead

at each time interval in the control (dcont), the number dead at

each time interval in the treatment (dexp), and the number of time

intervals that had valid data (m) *)



t=Table[input[[i,1]],{i,inc}];

dcont=Table[input[[i,2]],{i,inc}];

dexp=Table[input[[i,3]],{i,inc}];

m=inc;

pdeadc=dcont[[m]]/ncont;



(* Calculate the number of insects that die during each time

interval (kexp or kcont), the proportion that die during each time

interval (pexp or pcont), and correct for 0 dead. *)



kexp=Table[If[i==1,dexp[[1]],dexp[[i]]-dexp[[i-1]]],{i,m}];

pexpinit=kexp/nexp //N;

pexpcorr=Table[If[pexpinit[[i]]==0,0.5/(nexp+1),pexpinit[[i]]],{i,m}];

kcont=Table[If[i==1,dcont[[1]],dcont[[i]]-dcont[[i-1]]],{i,m}];

pcont=Table[kcont[[i]]/ncont,{i,m}];



(* Calculate the covariance matrix of p. *)



covp=Table[If[i==j,(pexpcorr[[i]]*(1-pexpcorr[[i]]))/

  (nexp-1),-(pexpcorr[[i]]*pexpcorr[[j]])/(nexp-1)],{i,m},{j,m}];



(* Apply Abbott's formula to correct for control mortality.  

  Calculate covariance matrix for corrected p's. *)



pexp2=Table[(pexpinit[[i]]-pcont[[i]])/(1-pcont[[i]]),{i,m}];

pexp=Table[If[pexp2[[i]]<=0,0.5/(nexp+1),pexp2[[i]]],{i,m}];



covpcorr=Table[If[i==j,((1/(1-pcont[[i]]))^2)*(((1-pexp[[i]])*

  pexp[[i]])/(nexp-1))+((1-pexp[[i]])/((1-pcont[[i]])^2)^2)*

  (((1-pcont[[i]])*pcont[[i]])/(ncont-1)),(1/(1-pcont[[i]]))*

  (1/(1-pcont[[j]]))*(-(pexp[[i]]*pexp[[j]])/(nexp-1))+

  (-(1-pexp[[i]])/((1-pcont[[i]])^2))*(-(1-pexp[[j]])/

  ((1-pcont[[j]])^2))*(-(pcont[[i]]*pcont[[j]])/(ncont-1))],

  {i,m},{j,m}];



(* Calculate the c and x matrices. *)



c=Table[If[j<=i,1,0],{i,m},{j,m}];

x=Table[If[j==1,1,t[[i]]],{i,m},{j,2}];



(* Calculate the cumulative proportion killed during each interval 

(pstarobs) and the covariance matrix of pstarobs. *)



pstarobs=c.pexp;

covpstar=c.covpcorr.Transpose[c];



(* If using probit transform, calculate the probits (z) of the

cumulative proportion killed during each interval, the probability

density of the probits (phiz), and the covariance matrix of the

probits (covz).  If using logits calculate the logits (z) and the

covariance matrix of the logits (covz).  If using CLL

transformation, calculate the transforms (z) and the covariance

matrix of the transforms. *)



Do[If[trany != \"probit\",Goto[other],y=Map[Quantile[ndist,#]&,pstarobs]];

z=Table[y[[i]],{i,m}];

phiz=(E^(-(z^2)/2))/Sqrt[2*Pi] //N;

covz=Table[If[i==j,covpstar[[i,i]]/phiz[[i]]^2,

  covpstar[[i,j]]/(phiz[[i]]*phiz[[j]])],{i,m},{j,m}];

Goto[endtrans];

Label[other];

If[trany != \"logit\",Goto[cll],z=Table[Log[pstarobs[[i]]/

  (1-pstarobs[[i]])],{i,m}]];

covz=Table[If[i==j,1/(nexp*pstarobs[[i]]*(1-pstarobs[[i]])),

  (1/(pstarobs[[i]]*(1-pstarobs[[i]])))*(1/(pstarobs[[j]]*

  (1-pstarobs[[j]])))*covpstar[[i,j]]],{i,m},{j,m}];

Goto[endtrans];

Label[cll];

z=Table[Log[-Log[1-pstarobs[[i]]]],{i,m}];

covz=Table[If[i==j,pstarobs[[i]]/(nexp*(1-pstarobs[[i]])*

  ((Log[1-pstarobs[[i]]])^2)),(1/(1-pstarobs[[i]]))*

  (1/-Log[1-pstarobs[[i]]])*(1/(1-pstarobs[[j]]))*

  (1/-Log[1-pstarobs[[j]]])*covpstar[[i,j]],],{i,m},{j,m}];

Label[endtrans];

,{i2,1}];



(* Calculate the (X'X)-1 (binit1) and X'Y (binit2) matrices.

Calculate the slope and intercept (the matrix newb) as binit1 times

binit2. *)



binit1=Inverse[Transpose[x].Inverse[covz].x];

binit2=(Transpose[x].Inverse[covz].z);

newb=binit1.binit2;

g1=ListPlot[Table[{t[[i]],z[[i]]},{i,m}],DisplayFunction

  ->Identity];



(* Begin iterations to converge on best values of slope and

intercept.  Calculate new probits (or logits or CLL

transforms) using slope and intercept (y).  Calculate new cumulative

proportion dead at each interval (pstar) from probits (or logits or

CLL transforms).  Calculate new proportion dying during each

interval (p) from pstar.  Calculate covariance matrix of p (covp)

and covariance matrix of pstar (covpstar).  Calculate new

probability density of probits (phiz), if using probit

transformation, and covariance matrix of the probits (or logits or

CLL transform) (covz).  Obtain new estimates of slope and

intercept (newb), and keep iterating until convergence criterion

met. *)



Do[

b=newb;

y=b[[1]]+(b[[2]]*t);

z=Table[y[[i]],{i,m}];

If[trany != \"probit\",Goto[other2],pstar=Map[CDF[ndist,#]&, z]];

p=Table[If[i==1,pstar[[1]],pstar[[i]]-pstar[[i-1]]],{i,m}];

covp=Table[If[i==j,(p[[i]]*(1-p[[i]]))/(nexp-1),

  -(p[[i]]*p[[j]])/(nexp-1)],{i,m},{j,m}];

phiz=(E^(-(z^2)/2))/Sqrt[2*Pi] //N;

covpstar=c.covp.Transpose[c];

covz=Table[If[i==j,covpstar[[i,i]]/phiz[[i]]^2,

  covpstar[[i,j]]/(phiz[[i]]*phiz[[j]])],{i,m},{j,m}];

Goto[endtrans2];

Label[other2];

If[trany != \"logit\",Goto[cll2],pstar=Exp[z]/(1+Exp[z])];

p=Table[If[i==1,pstar[[1]],pstar[[i]]-pstar[[i-1]]],{i,m}];

covp=Table[If[i==j,(p[[i]]*(1-p[[i]]))/(nexp-1),

  -(p[[i]]*p[[j]])/(nexp-1)],{i,m},{j,m}];

covpstar=c.covp.Transpose[c];

covz=Table[If[i==j,1/(nexp*pstar[[i]]*(1-pstar[[i]])),

  (1/(pstar[[i]]*(1-pstar[[i]])))*(1/(pstar[[j]]*

  (1-pstar[[j]])))*covpstar[[i,j]]],{i,m},{j,m}];

Goto[endtrans2];

Label[cll2];

pstar=1-Exp[-Exp[z]];

p=Table[If[i==1,pstar[[1]],pstar[[i]]-pstar[[i-1]]],{i,m}];

covp=Table[If[i==j,(p[[i]]*(1-p[[i]]))/(nexp-1),

  -(p[[i]]*p[[j]])/(nexp-1)],{i,m},{j,m}];

covpstar=c.covp.Transpose[c];

covz=Table[If[i==j,pstar[[i]]/(nexp*(1-pstar[[i]])*

  ((Log[1-pstar[[i]]])^2)),(1/(1-pstar[[i]]))*

  (1/-Log[1-pstar[[i]]])*(1/(1-pstar[[j]]))*

  (1/-Log[1-pstar[[j]]])*covpstar[[i,j]],],{i,m},{j,m}];

Label[endtrans2];

binit1=Inverse[Transpose[x].Inverse[covz].x];

binit2=(Transpose[x].Inverse[covz].z);

newb=binit1.binit2;

If[Abs[newb[[1]]-b[[1]]]<.00000001 && 

  Abs[newb[[2]]-b[[2]]]<.00000001,Break[]];,{i,50}];



(* Calculate the standard errors of the slope (seslope) and intercept

(seint) from the covariance matrix of b (binit1). *)



seint=Sqrt[binit1[[1,1]]];

seslope=Sqrt[binit1[[2,2]]];



(* Write output about slope and intercept to file. *)



WriteString[out,\"Intercept = \"];

Write[out,newb[[1]]];

WriteString[out,\"\\nSE Intercept = \"];

Write[out,seint];

WriteString[out,\"\\nSlope = \"];

Write[out,newb[[2]]];

WriteString[out,\"\\nSE Slope = \"];

Write[out,seslope];

WriteString[out,\"\\nCovariance of Slope & Intercept = \"];

Write[out,binit1[[1,2]]];

WriteString[out,\" \\n\"];



WriteString[out,

\"The following five lines can be used to enter slopes and intercepts 

into the program RELPOT for calculating relative potency or into the 

program SLOPE for calculating equality of slopes and intercepts\"];

WriteString[out,\"\\n\\na = \"];

Write[out,newb[[1]]];

WriteString[out,\"sea = \"];

Write[out,seint];

WriteString[out,\"b = \"];

Write[out,newb[[2]]];

WriteString[out,\"seb = \"];

Write[out,seslope];

WriteString[out,\"cov = \"];

Write[out,binit1[[1,2]]];

WriteString[out,\" \\n\"];





(* Calculate Chi-square to test for goodness of fit of regression

line to observed data.  Probit, logit, or CLL transform

(zpred) is calculated from final values of slope and intercept.  A

new p (ppred) is calculated from the probit (or logit or CLL

transform).  The cumulative number dead at each time interval (dpred) is

calculated from the p and ncorr (the number of insects tested minus

the number that died in the control).  The noncumulative number

dying during each time interval (kpred) is calculated from dpred.

The observed number dead at each time interval (kpred) is calculated

from the cumulative proportion dead at each time interval (pexp)

times the initial number of treated insects (ncorr).  The observed

minus expected number of insects dead is calculated as dobs - dexp.

The covariance matrix for the observed minus expected number dead

(covobsme) is calculated as ncorr squared times the covariance

matrix of p (covp).  The Chi-square value (chisq) is then

calculated, and the probability of a larger Chi-square (prchisq)

with n-2 degrees of freedom (dfchisq) is determined. *)



zpred=newb[[1]]+(newb[[2]]*t);

Do[

If[trany != \"probit\",Goto[other3],ppred=Map[CDF[ndist,#]&,zpred]];

Goto[endtrans3];

Label[other3];

If[trany != \"logit\",Goto[cll3],ppred=Exp[zpred]/(1+Exp[zpred])];

Goto[endtrans3];

Label[cll3];

ppred=1-Exp[-Exp[zpred]];

Label[endtrans3];

,{i3,1}];

ncorr=nexp*(1-pdeadc);

dpred=ppred*ncorr;

kpred=Table[If[i==1,dpred[[1]],dpred[[i]]-dpred[[i-1]]],{i,m}];

kobs=pexp*ncorr;

obsmexp=kobs-kpred;

covobsme=(ncorr^2)*covp;

chisq=obsmexp.Inverse[covobsme].obsmexp;

chidist=ChiSquareDistribution[m-2];

prchisq=2*PDF[chidist,chisq];

dfchisq=m-2;



(* If the Chi-square is significant, correct all variances used to

calculate confidence limits on LT levels by multiplying each variance

by h. *)



h=If[prchisq<.05,chisq/(m-1),1];



(* Chi-square information is output to file. *)



WriteString[out,\"Obs - Exp number of insects dying during each time

interval = \\n\"];

Write[out,obsmexp];

WriteString[out,\"\\nChiSquare = \"];

Write[out,chisq];

WriteString[out,\"\\nProbability of a greater ChiSquare = \"];

Write[out,FortranForm[prchisq]];

WriteString[out,\"\\nDegrees of freedom for ChiSquare = \"];

Write[out,dfchisq];

WriteString[out,\" \\n\"];

If[prchisq>0.05,Null,WriteString[out,\"Heterogeneity factor = \"];

  Write[out,h]];

WriteString[out,\" \\n\"];

g2=Plot[Table[newb[[1]]+(newb[[2]]*i)],{i,input[[1,1]],t[[m]]},

  DisplayFunction->Identity];



If[graphs==\"no\",Goto[nographs]];

Do[

If[trany != \"probit\",Goto[other4]];

If[tranx == \"yes\",Goto[logyes1],Show[g1,g2, DisplayFunction -> 

  (Print[#]; &),Frame->True,PlotLabel-> 

  \"Observed and Calculated Probits vs Time\",DisplayFunction->

  $DisplayFunction,AxesOrigin->{0,0}]];

Goto[endtrans4];

Label[logyes1];

Show[g1,g2, DisplayFunction -> (Print[#]; &),Frame->True,PlotLabel->

  \"Observed and Calculated Probits vs Log Time\",DisplayFunction->

  $DisplayFunction,AxesOrigin->{0,0}];

Goto[endtrans4];

Label[other4];

If[trany != \"logit\",Goto[cll4]];

If[tranx == \"yes\",Goto[logyes2],Show[g1,g2, DisplayFunction -> 

  (Print[#]; &),Frame->True,PlotLabel->

  \"Observed and Calculated Logits vs Time\",DisplayFunction->

  $DisplayFunction,AxesOrigin->{0,0}]];

Goto[endtrans4];

Label[logyes2];

Show[g1,g2, DisplayFunction -> (Print[#]; &),Frame->True,PlotLabel->

  \"Observed and Calculated Logits vs Log Time\",DisplayFunction->

  $DisplayFunction,AxesOrigin->{0,0}];

Goto[endtrans4];

Label[cll4];

If[tranx == \"yes\",Goto[logyes3],Show[g1,g2, DisplayFunction -> 

  (Print[#]; &),Frame->True,PlotLabel->

  \"Observed and Calculated CLL transforms vs Time\", DisplayFunction->

  $DisplayFunction,AxesOrigin->{0,0}]];

Goto[endtrans4];

Label[logyes3];

Show[g1,g2, DisplayFunction -> (Print[#]; &),Frame->True,PlotLabel->

\"Observed and Calculated CLL transforms vs Log Time\", DisplayFunction->

  $DisplayFunction,AxesOrigin->{0,0}];

Label[endtrans4];

,{i4,1}];

Label[nographs];



(* Use Fieller's theorem to calculate confidence limits on LT levels.

 *)



tdist=StudentTDistribution[m-2];

t2=If[prchisq>.05,1.96,Quantile[tdist,.975]];

If[prchisq>.05,WriteString[out,\"z for calculating confidence limits =

\"],WriteString[out,\"t for calculating confidence limits = \"]];

Write[out,t2];

WriteString[out,\" \\n\"];



a=newb[[1]];

b=newb[[2]];

v11=binit1[[1,1]]*h;

v22=binit1[[2,2]]*h;

v12=binit1[[1,2]]*h;



outprob1=Table[i,{i,.1,.9,.1}];

outprob2=Table[i,{i,.91,.99,.01}];

outprob3=Table[i,{i,.991,.999,.001}];



Do[

If[trany != \"probit\",Goto[other5],lt1=Map[Quantile[ndist,#]&,

  outprob1]];

lt2=Map[Quantile[ndist,#]&,outprob2];

lt3=Map[Quantile[ndist,#]&,outprob3];

Goto[endtrans5];

Label[other5];

If[trany != \"logit\",Goto[cll5],lt1=Table[Log[outprob1[[i]]/

  (1-outprob1[[i]])],{i,9}]];

lt2=Table[Log[outprob2[[i]]/(1-outprob2[[i]])],{i,9}];

lt3=Table[Log[outprob3[[i]]/(1-outprob3[[i]])],{i,9}];

Goto[endtrans5];

Label[cll5];

lt1=Table[Log[-Log[1-outprob1[[i]]]],{i,9}];

lt2=Table[Log[-Log[1-outprob2[[i]]]],{i,9}];

lt3=Table[Log[-Log[1-outprob3[[i]]]],{i,9}];

Label[endtrans5];

,{i5,1}];



m1=(lt1-a)/b;

m2=(lt2-a)/b;

m3=(lt3-a)/b;

x11=m1+(v12/v22);

x12=m2+(v12/v22);

x13=m3+(v12/v22);

x21=v11+2*m1*v12+(m1^2)*v22;

x22=v11+2*m2*v12+(m2^2)*v22;

x23=v11+2*m3*v12+(m3^2)*v22;

x3=v11-((v12^2)/v22);

g=(t2^2)*v22/(b^2);

y11=m1+(g*x11/(1-g));

y12=m2+(g*x12/(1-g));

y13=m3+(g*x13/(1-g));

y21=t2*Sqrt[x21-(g*x3)]/(b*(1-g));

y22=t2*Sqrt[x22-(g*x3)]/(b*(1-g));

y23=t2*Sqrt[x23-(g*x3)]/(b*(1-g));

l1=y11-y21;

l2=y12-y22;

l3=y13-y23;

u1=y11+y21;

u2=y12+y22;

u3=y13+y23;

Do[If[g*x3<x21[[i]] && (1-g)>0,Goto[Skip1]];

l1[[i]]=\"Undefined\";

u1[[i]]=\"Undefined\";

Label[Skip1],{i,9}];

Do[If[g*x3<x22[[i]] && (1-g)>0,Goto[Skip2]];

l2[[i]]=\"Undefined\";

u2[[i]]=\"Undefined\";

Label[Skip2],{i,9}];

Do[If[g*x3<x23[[i]] && (1-g)>0,Goto[Skip3]];

l3[[i]]=\"Undefined\";

u3[[i]]=\"Undefined\";

Label[Skip3],{i,9}];

Do[If[m1[[i]]<0,m1[[i]]=0],{i,9}];

Do[If[m2[[i]]<0,m2[[i]]=0],{i,9}];

Do[If[m3[[i]]<0,m3[[i]]=0],{i,9}];

Do[If[l1[[i]]<0,l1[[i]]=0],{i,9}];

Do[If[l2[[i]]<0,l2[[i]]=0],{i,9}];

Do[If[l3[[i]]<0,l3[[i]]=0],{i,9}];

Do[If[u1[[i]]<0,u1[[i]]=0],{i,9}];

Do[If[u2[[i]]<0,u2[[i]]=0],{i,9}];

Do[If[u3[[i]]<0,u3[[i]]=0],{i,9}];



(* Output results of Fieller's theorem calculations to file. *)



If[tranx==\"yes\",WriteString[out,\"Log10 transform was used; following

data output as logarithms\\n\"]];

WriteString[out,\"\\n   Prop.      Lower      Time       Upper\\n\"];

WriteString[out,\"   kill       limit                 limit\\n\"];

WriteString[out,\" \\n\"];

Do[Write[out,PaddedForm[SequenceForm[outprob1[[i]],\"  \",l1[[i]],

  \"  \",m1[[i]],\"  \",u1[[i]]],{7,4}]],{i,9}];

Do[Write[out,PaddedForm[SequenceForm[outprob2[[i]],\"  \",l2[[i]],

  \"  \",m2[[i]],\"  \",u2[[i]]],{7,4}]],{i,9}];

Do[Write[out,PaddedForm[SequenceForm[outprob3[[i]],\"  \",l3[[i]],

  \"  \",m3[[i]],\"  \",u3[[i]]],{7,4}]],{i,9}];



(* If using log transform, write out LT values in real time. *)



Do[

If[tranx==\"no\",Goto[end]];

outl1=Table[If[l1[[i]]==\"Undefined\",\"Undefined\",10^l1[[i]]],{i,9}];

outl2=Table[If[l2[[i]]==\"Undefined\",\"Undefined\",10^l2[[i]]],{i,9}];

outl3=Table[If[l3[[i]]==\"Undefined\",\"Undefined\",10^l3[[i]]],{i,9}];

outm1=Table[10^m1[[i]],{i,9}];

outm2=Table[10^m2[[i]],{i,9}];

outm3=Table[10^m3[[i]],{i,9}];

outu1=Table[If[u1[[i]]==\"Undefined\",\"Undefined\",10^u1[[i]]],{i,9}];

outu2=Table[If[u2[[i]]==\"Undefined\",\"Undefined\",10^u2[[i]]],{i,9}];

outu3=Table[If[u3[[i]]==\"Undefined\",\"Undefined\",10^u3[[i]]],{i,9}];



WriteString[out,\"\\n\\nOutput converted to original time units\\n\\n\"];

WriteString[out,\"   Prop.      Lower      Time       Upper\\n\"];

WriteString[out,\"   kill       limit                 limit\\n\\n\"];

Do[Write[out,PaddedForm[SequenceForm[outprob1[[i2]],\"  \",outl1[[i2]]

  ,\"  \",outm1[[i2]],\"  \",outu1[[i2]]],{7,4}]],{i2,9}];

Do[Write[out,PaddedForm[SequenceForm[outprob2[[i2]],\"  \",outl2[[i2]]

  ,\"  \",outm2[[i2]],\"  \",outu2[[i2]]],{7,4}]],{i2,9}];

Do[Write[out,PaddedForm[SequenceForm[outprob3[[i2]],\"  \",outl3[[i2]]

  ,\"  \",outm3[[i2]],\"  \",outu3[[i2]]],{7,4}]],{i2,9}];

Label[end];

,{i3,1}];



loop=loop+1;

If[loop<7,Goto[loopstart]];

,{iloop,1}];



Close[out]\
\>", "Input"]
},
WindowSize->{707, 787},
WindowMargins->{{24, Automatic}, {Automatic, 30}},
FrontEndVersion->"9.0 for Microsoft Windows (64-bit) (January 25, 2013)",
StyleDefinitions->"Default.nb"
]
(* End of Notebook Content *)

(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[557, 20, 19465, 1096, 18629, "Input"]
}
]
*)

(* End of internal cache information *)

