Symbolic and Monte Carlo
number keying error analysis

Harold Thimbleby
harold@thimbleby.net
Swansea University
United Kingdom

This file analyzes number keying errors, using both a symbolic approach and a Monte Carlo method. See the web site at http://www.cs.swan.ac.uk/~csharold/som for more details.

Start here…

To use this Mathematica notebook for the first time: do Evaluation>Evaluate Notebook (that is, in the Evaluation menu, select Evaluate Notebook). Everything should work, and you will be asked whether you want to try running some slow calculations or to skip them; I recommend skipping them!

You may wish to do Cell>Delete All Output if you are not familiar with Mathematica, as this will remove any possibly confusing generated results left over from a previous run.

Once you have run everything, the data can be saved in a file, and you can recover it by executing the Get command, which you can find below. If you use Get in another session to restore data, don't forget that you must also do all initialization again, else the rest of the code won't work properly! You can get Mathematica to do initialisation only by selecting Evaluation>Evaluate Initialization Cells, but you will still need to go through this notebook and evaluate cells to generate data, analyze it and visualize. To get familiar with what's going on, it's easier to Evaluation>Evaluate Notebook and do everything!

Note that you can combine various options of error analysis together (e.g., how number errors are handled) and choose what to visualize later. The data structures will keep track of what you did.

If you are lucky, you will get lots of data and graph visualizations. Have fun.

The following code assumes you have either selected every cell or you are going to evaluate every cell yourself. It explains what fastRun does and initializes it. It'll be run automatically if you Evaluate Notebook, as above.

fastRun=True;(* You may want to make this False next time you run this notebook! *)
If[fastRun,fastRun=ChoiceDialog["If this is your first run, you probably want a quick overview of illustrative results, rather than detailed results that take a lot longer to wait for.\n\nIf you select \"quick results\", by clicking the option below, you will quickly generate a small dataset, and you will be given a choice of symbolic or Monte Carlo visualizations (you can run the visualizations again on the same data without regenerating it; you can also add further data as you experiment).\n\nTo see results, you will need to scroll down the notebook to where they are displayed (and probably wait for Mathematica to catch up with you).",{"I want quick results"->True,"I want to wait for more detailed results"->False}]];

"Last run on "<>DateString["Date"]<>" at "<>DateString["Time"]<>" using Mathematica "<>$Version<>"\n\n"<>If[fastRun,"Quick results","Detailed results"]<>" selected."

How it works — overview

It's simplest to consider the Monte Carlo analysis first. Think of a number, then try keying it with certain probabilities of making keying errors. Treat that possibly erroneous key sequence like a typical gadget, and get a numeric value (which would be what a user would get using a typical device). The ratio of this value to the intended value will be r; thus there is an out by r error. Of course r may be 1, in which case there is no numerical error. We then classify the keyed number as syntactically valid (good) or invalid (bad). With a bit of tinkering, this allows us to plot graphs of the probability of an out by 10 error for various underlying probabilities of keying errors.

More precisely, for each number, the Monte Carlo analysis tries, for a range of probabilities (e=0.01 to 1.00), keying the number. This gives for each out by r error a vector of numbers, the frequencies those out by r errors occur in the simulation for each e. To be any use later in the analysis, these vectors are converted to polynomials in e.

Thinking of a random number is done as follows. Generate four decimal digits (from a uniform distribution), and convert to the form XX.XX. Remove leading and trailing zeros. This is now a number between 0.01 and 99.99 and in a valid form.

This Mathematica notebook also provides symbolic analyses. These are slower, and typically you won't want to wait around for a symbolic analysis of 0.01 to 99.99. The advantage is that the analysis is exact; here is how it is done.

Instead of choosing a single random number, for every valid number in a given format (e.g., 0.1 to 9.9), consider every way of keying it (rather than just one way, chosen probabalistically). For each way of keying it, work out the probability that the user would have keyed it that way. We know the underlying probability rate, e, and this calculation is purely symbolic. As with the Monte Carlo approach, it ends up with polynomials in e, but this time they are exact.

Assumptions

It could take forever if we model the user keying anything. Instead, when analyzing how the user keys a number, we only consider the user keying the correct number of keystrokes. We allow for the user prematurely terminating the number, however (because we treat termination as a keystroke with a certain probability), but we do not analyze if the user keys too much. In fact, it's clear that if a keyed number is out by 10, keying more digits isn't going to stop it being out by ten. However, this assumption in our models means that we will miss some out by 10 errors.

We do not define a number with too many keys as an invalid number. We are only looking at numbers 0.01 to 99.99, so this is not an issue.

We do not require the use of commas in numbers ≥1,000 for similar reasons.

We treat a single zero (0) as an error as a drug dose of zero would be incorrect, but this isn't an explicit ISMP rule.

Perhaps, ideally, the analysis could be more realistic, but it would just make it more complex without any obvious argument as to the benefits of better understanding what's going on.

The Monte Carlo analysis has a hard-wired number of runs. It is set very small for a first run (100?), and larger otherwise (10,000?), but you can change it easily by editing the code.

Defining probabilities

Probability of hitting [0123456789.] in error is p; probability of hitting $ (representing the termination key, perhaps called ENTER) in error is q. We assume all digits/dot equiprobable.

We define k=p/q (k would be 2 for example if the user is twice as likely to press the wrong numeric key than end the number prematurely), and then for an overall error probability e=p+q, e=[0,1], we get:

Clear[e,p,q,k];
Solve[e==p+q&&k==p/q,{q,p}]/.(a_->b_):>(a=b);
Grid[{{Defer[p],"=",p},{},{Defer[q],"=",q}}]//TraditionalForm

For my iPhone calculator, there are 11 digit keys and 4 ways to end a number, so we estimate k=7/4; on the Abbott aimplus pump, there are 11 digit keys and only 1 key ends a number, called [YES/ENTER] on the Abbott aimplus, so we estimate k=(11-1)/1:

k=10;

Auxilliary functions etc

The obvious way to represent a keyed number is as a list of its keystrokes. A keystroke is coded as a number, with 0–9 being the digit keys 0–9, and decimal point being coded as 10, and the end of a number (which we call dollar) being coded as 11.

It's convenient to define names for the decimal and terminating codes, and a convenience function to convert one of these number lists into a printable string.

dot=10;
dollar=11;
asString[list_]:=StringJoin@@(Which[#==dot,".",#==dollar,"$",True,ToString[#]]&/@list);

For example,

asString@{1,2,3,dot,9,0,dollar}//FullForm

We need to know the value of a number and whether it is a syntax error. We define several ways of reading numbers…

Note that when you generate data, the dataset records which version of value[] was used, so you can come back here, change the error handling assumptions and generate some more data if you wish. (However, as currently configured, although you can collect as much data as you wish, you can't visualize different sorts of data together.)

All versions of the code assume the only possible keys are 0123456789, dot or dollar.

Version 1. Ignore errors (common on calculators)

Here, value[] keeps on trying to read the number regardless of errors; so 1.2.3 is returned as 1.23.

valueIgnore[list_]:=(*parse number; stops at end, or non-digit after first decimal point*)
Module[{c,n=0,len=Length[list],d},
For[i=1,i≤len,i++,
c=list[[i]];
If[0≤c&&c≤9,n=n*10+c,Break[]]
];
If[c==dot,
For[d=10;i++,i≤len,i++,
c=list[[i]];
If[0≤c&&c≤9,n=n+c/d;d=d*10,
If[c==dollar,Break[]]
]
]
];
N[n]
]

Version 2. Stop at first error (common on infusion pumps)

Here, value[] stops at the first error; so 1.2.3 is returned as 1.2.

valueStop[list_]:=
Module[{c,n=0,len=Length[list],d},
For[i=1,i≤len,i++,
c=list[[i]];
If[0≤c&&c≤9,n=n*10+c,Break[]]
];
If[c==dot,
For[d=10;i++,i≤len,i++,
c=list[[i]];
If[0≤c&&c≤9,n=n+c/d;d=d*10,
Break[]
]
]
];
N[n]
]

Version 3. Decimal point resets decimal part to zero (like Graseby 3400 infusion pump)

Here, value[] decimal point clears the decimal part; so 1.2.3 is returned as 1.3.

valueGraseby[list_]:=
Module[{c,n=0,f=0,len=Length[list],d},
For[i=1,i≤len,i++,
c=list[[i]];
If[0≤c&&c≤9,n=n*10+c,Break[]]
];
If[c==dot,
For[d=10;i++,i≤len,i++,
c=list[[i]];
If[0≤c&&c≤9,f=f+c/d;d=d*10,
If[c==dollar,Break[],f=0;d=10]
]
]
];
N[n+f]
]

Version 4. Use last decimal point (logical idea, but haven't yet seen in used)

Here, value[] takes the relevant decimal point as the last one; so 1.2.3 is returned as 12.3. The easiest way to do this would be to delete the earlier decimal points, then call one of the other value[] functions. But this code is called frequently, and it's efficiency is important, so we've written it in a less-clear way that should be faster.

valueLater[list_]:=
Module[{c,n=0,nn,f=0,len=Length[list],d},
For[i=1,i≤len,i++,
c=list[[i]];
If[0≤c&&c≤9,
n=n*10+c,
Break[];
]
];
If[c==dot,
nn=n;
d=10;
For[i=i+1,i≤len,i++,
c=list[[i]];
If[0≤c&&c≤9,
nn=nn*10+c;
f=f+c/d;
d=d*10,
If[c==dot,
n=nn;
f=0;
d=10,
Break[]
]
];
]
];
N[n+f]
]

{value,valueType}=ChoiceDialog["What sort of number handling do you want analyzed?",{"Use first decimal point\n1.2.3 is treated as 1.23"->{valueIgnore,"use first decimal point\n1.2.3->1.23"},"Stop at errors\n1.2.3 is treated as 1.2"->{valueStop,"stop at errors: 1.2.3->1.2"},"Decimal point resets decimal part\n1.2.3 is treated as 1.3"->{valueGraseby,"decimal points reset: 1.2.3->1.3"},"Use last decimal point\n1.2.3 is treated as 12.3"->{valueLater,"use last decimal point: 1.2.3->12.3"}}];

caption[]="Analysis type: "<>valueType;

Module[{example={1,2,dot,2,6,dot,3,4,dollar,9}},"You selected \""<>valueType<>"\" for handling number errors.\nHere's an example how bad numbers are handled with your choice: "<>asString[example]<>" -> "<>ToString@value[example]
]

Define error function

error[] returns true if the number (perhaps truncated by a $) is invalid by ISMP rules.

ISMP rules are: if an integer (no decimal part) there should be no leading decimal; if a decimal, then no leading or trailing zeros, with the exception of 0.XXXX; and there should be no naked decimals, that is the decimal can't be the first or last character. To which I add at most one decimal point, which ISMP don't mention!

The code is written in this way to make it a bit more efficient that a more declartive style — the code is used in lots of inner loops, so its speed is important.

error[v_]:=
Module[{n,dots=0,max=Length[v],first,last},
For[n=1,n≤max,n++,
If[v[[n]]==dot,dots++];
If[v[[n]]==dollar,max=n-1;Break[]]
];
If[max==0,Return[True]];(* an empty list is an error *)
If[dots>1,Return[True]];(* multiple decimal points is an error *)
first=v[[1]];
If[dots==0,Return[first==0]];(* if an integer, only problem is leading zero (includes equal to zero) *)

(* remaining tests are for numbers with exactly 1 decimal point *)
If[first==dot,Return[True]];(* must not have leading decimal point *)
last=v[[max]];
(* if has 1 decimal point and it's not first or last, then there are at least 3 characters *)
last==dot||last==0||(first==0&&v[[2]]≠dot)
]

For example, here's a table of various numbers:

Grid@Map[Map[Function[s,Style[s,If[error[#],Red,Darker[Green]]]],{asString[#],If[error[#]," invalid","ok"],"=",value[#]}]&,{{1,dot,2},{1,dot,2,dot,3},{1,dot,2,dollar,9,dot},{dot,2},{0,dot,1},{dollar},{6,0,dot,0},{6,0,dot,1,0},{0,6,dot,1},{6,0,dot,1},{1,2,0},{0,1,2},{1,dot,dollar},{0}}]

Define bins

The dose error is the ratio of the intended and actual keyed doses. Call the ratio t; we need to track the probability of each event t, which reprsents an out by t error. To program that we keep an array b with scale bins, such that b[i] will accummulate the probability of out by t events where Round[scale*t]=i. Since Mathematica arrays have upper bounds, any t≥huge/scale is mapped to huge. For all our graph plotting, we're only really interested in plots to t=10, so huge is set to 12scale. (For our purposes we just need a large enough value of scale to make the graphs look smooth.)

scale=100;
huge=12scale;
makeArray[]:=makeArray[0];
makeArray[n_]:=Table[n,{huge}];

With a scale of 100, the bins will be 1/scale=0.01 wide.

Having got the bins, we may want to find values in them! Given lists of unsorted {x,y} data, interpolate to find the point closest to a given x or y. Naff algorithm, but it works reliably with few assumptions.

xintercept[n_,list_]:=Module[{lo,hi},
lo=Last@Sort[Cases[list,{_,_?(#≤n&)}],#1[[2]]<#2[[2]]&];
hi=First@Sort[Cases[list,{_,_?(#≥n&)}],#1[[2]]<#2[[2]]&];
If[lo[[2]]==n,{lo[[1]],n},{lo[[1]]+(hi[[1]]-lo[[1]])(n-lo[[2]])/(hi[[2]]-lo[[2]]),n}]
]

yintercept[n_,list_]:=Module[{lo,hi},
lo=Last@Sort[Cases[list,{_?(#≤n&),_}],#1[[1]]<#2[[1]]&];
hi=First@Sort[Cases[list,{_?(#≥n&),_}],#1[[1]]<#2[[1]]&];
If[lo[[1]]==n,{n,lo[[2]]},{n,lo[[2]]+(hi[[2]]-lo[[2]])(n-lo[[1]])/(hi[[1]]-lo[[1]])}]
]

Example:

xintercept[4,{{1,2},{2,3},{4,5},{6,7}}]

Generate data

The code works by generating all numbers of a particular correct syntactic form (say, 10 to 99), and for each such number (say 23) generating all possible ways of keying it. This is done by generating all base 12 numbers and converting them to lists — these are then all possible key sequences, usually called dj in the code. For each key sequence we work out its value and (given p and q) the probability that the user would have keyed it.

Each syntactic form keeps to arrays of bins: good and bad (or g and b in the code). It would be nice to track out by t errors, but that would be slow; instead we track out by "exactly" t, to within the precision given by the scale of the bins.

Good tracks the out by 'exactly' t errors that are not syntax errors, and bad tracks the out by t errors that are syntax errors. All the functions below then returns the pair of arrays (along with some other descriptive stuff) an appends to the vector data.

Make some constants, so we don't keep repeating these sums in inner loops:

ok=1-p-q;
poverten=p/10;

SetAttributes[update,HoldAll];
update[dj_,b_,g_,val_,prob_]:=
(* 
dj is what was keyed,
b and g are the bad and good bin arrays,
val is the value of the intended number,
prob is the probability of keying dj 
*)
Module[
{t=value[dj]},(* t is the device's idea of the value of the keyed number *)

(* set t to the scaled out by t ratio so it can be used directly as an index into the bins; if the scaled t is out of range, make it huge *)
t=If[t==0||(t=scale*Max[val/t,t/val])≥huge,huge,Round[t]];

(* if the keyed number is invalid, add to bad bin, else add to good bin *)
If[error[dj],b[[t]]+=prob,g[[t]]+=prob]
];

Clear data for a new run

Clear data to initialize for a new run…

data={};
dataset=0;(* give each data set a unique number - helps if you rerun things *)

Symbolic: 1 digit integer, 1-9

Having done all that, the easiest code to explain in detail is to analyze all possible ways of keying a single digit. Please see all the comments below; all further symbolic approaches use the same method.

AppendTo[data,
Module[{g,b,i,j,dj,prob,val,entries=0},
g=makeArray[];(* generate empty bins *)
b=makeArray[];
Print@ProgressIndicator[Dynamic[i/10]];
For[i=1,i<10,i++,(* generate correct numbers 1 to 9 *)
val=i;(* the value of a single digit is easy: it's itself *)
entries++;
For[j=0,j<12,j++,(* all possible ways j of keying it *)
(* calculate the probability of keying this version *)
prob=
If[
j==dollar,q,(* dollar is always unlike dj[k], as there is no $ in the generated number *)
If[j≠i,poverten,(* poverten is a constant, p/10; there are 10 ways of being *unlike* the key it should be, 0-9 or dot (here, the right key is j) apart from $ *)
ok (* the digit keyed is correct, with probability ok *)
]
];
dj={j};(* for consistency, convert j to {j} which is the format we use everywhere else below *)
update[dj,b,g,val,prob](* update the good and bad bins *)
]
];
{"good"->g,"bad"->b,"description"->"1 digit with no decimal point, 1-9","size"->entries,"dataset"->++dataset,"type"->"Partial","analysis"->caption[]}
]
];

Symbolic: 2 digit integer, 10-99

AppendTo[data,
Module[{g,b,i,di,j,dj,prob,k,val,entries=0},
g=makeArray[];
b=makeArray[];
Print@ProgressIndicator[Dynamic[i/90]];
For[i=10,i<100,i++,
di=IntegerDigits[i];(* what user should key *)
val=i;
entries++;
For[j=0,j<12^2,j++,
dj=PadLeft[IntegerDigits[j,12],2]; (* what they do key *)
(* what is prob of keying dj? *)
prob=1;
For[k=1,k≤2,k++,
prob*=
If[
dj[[k]]==dollar,q,(* dollar is always unlike dj[k], as there is no $ in the generated number *)
If[dj[[k]]≠di[[k]],poverten,
ok
]
]
];
update[dj,b,g,val,prob]
]
];
{"good"->g,"bad"->b,"description"->"2 digit integer, 10-99","size"->entries,"dataset"->++dataset,"type"->"Partial","analysis"->caption[]}
]
];

Symbolic: 2 digits with a decimal point, 0.1-9.9 (obviously excluding 1.0, 2.0 which are invalid etc)

Unlike the previous examples, converting i to the numeric value isn't so trivial; we use value[] as defined above.

If[!fastRun,
AppendTo[data,
Module[{g,b,i,di,j,dj,prob,k,val,entries=0},
g=makeArray[];
b=makeArray[];
Print@ProgressIndicator[Dynamic[i/90]];
For[i=1,i<10^2,i++,
If[Mod[i,10]==0,Continue[]];(* easy way to ensure doesn't end in zero *)
di=PadLeft[IntegerDigits[i],2];(* what user should key *)
di={di[[1]],dot,di[[2]]};
val=value[di];
entries++;
For[j=0,j<12^3,j++,
dj=PadLeft[IntegerDigits[j,12],3]; (* what they do key *)
(* what is prob of keying dj? *)
prob=1;
For[k=1,k≤3,k++,
prob*=
If[
dj[[k]]==dollar,q,(* dollar is always unlike dj[k], as there is no $ in the generated number *)
If[dj[[k]]≠di[[k]],poverten,
ok
]
]
];
update[dj,b,g,val,prob]
]
];
{"good"->g,"bad"->b,"description"->"2 digits with a decimal point, 0.1-9.9","size"->entries,"dataset"->++dataset,"type"->"Partial","analysis"->caption[]}
]
]
];

Symbolic: 3 digits with a decimal point after second digit, 10.1-99.9 (obviously excluding 10-99 integers)

If[!fastRun,
AppendTo[data,
Module[{g,b,i,di,j,dj,prob,k,t,val,entries=0},
g=makeArray[];
b=makeArray[];
Print@ProgressIndicator[Dynamic[i/1000]];
For[i=101,i<1000,i++,
di=IntegerDigits[i];(* what user should key *)
If[di[[3]]==0,Continue[]];
di={di[[1]],di[[2]],dot,di[[3]]};
val=value[di];
entries++;
For[j=0,j<12^3,j++,
dj=PadLeft[IntegerDigits[j,12],3]; (* what they do key *)
(* what is prob of keying dj? *)
prob=1;
For[k=1,k≤3,k++,
prob*=
If[
dj[[k]]==dollar,q,(* dollar is always unlike dj[k], as there is no $ in the generated number *)
If[dj[[k]]≠di[[k]],poverten,
ok
]
]
];
update[dj,b,g,val,prob]
]
];
{"good"->g,"bad"->b,"description"->"3 digits with a decimal point after second digit, 10.1-99.9","size"->entries,"dataset"->++dataset,"type"->"Partial","analysis"->caption[]}
]
]
];

Symbolic: 1 to 4 digits with or without a decimal point, 0.01-99.99 (takes ages!)

This code, below, is optionally disabled because it is soooo slow: it's going to take much longer than all of the other bits of code put together. You can thus safely execute everything in this file, and it will work in a reasonable time (like, 10 minutes) — provided you answer the ChoiceDialog box question wisely!

If[!fastRun&&ChoiceDialog["Handling 0.01-99.99 symbolically is very, very slow.\n\nDo you really want to do it (and perhaps wait days)?",{"No. Skip this step."->False,"Yes. Wait ages."->True}],
AppendTo[data,
Module[{g,b,i,di,j,dj,k,prob,val,entries=0,cases,keyPresses},
g=makeArray[];
b=makeArray[];
Print@ProgressIndicator[Dynamic[i/10000]];
For[i=1,i<10000,i++,
di=PadLeft[IntegerDigits[i],4];
(* remove trailing zeros after decimal point; and remove decimal point if two trailing zeros *)
If[di[[4]]≠0,
di={di[[1]],di[[2]],dot,di[[3]],di[[4]]},
If[di[[3]]≠0,
di={di[[1]],di[[2]],dot,di[[3]]},
If[di[[1]]≠0,
di={di[[1]],di[[2]]},
di={di[[2]]}
]
]
];
If[di[[1]]==0,di=Rest[di]];(* delete any leading zero *)
(* di is now a list of what user should key *)
keyPresses=Length[di];
cases=12^keyPresses;
val=value[di];
entries++;
For[j=0,j<cases,j++,
dj=PadLeft[IntegerDigits[j,12],keyPresses]; (* what is actually keyed *)
(* what is prob of keying dj? *)
prob=1;
For[k=1,k≤keyPresses,k++,
prob*=
If[
dj[[k]]==dollar,q,
If[dj[[k]]≠di[[k]],poverten,
ok
]
]
];
update[dj,b,g,val,prob]
]
];
{"good"->g,"bad"->b,"description"->"upto 4 digits with or without a decimal point, 0.01-99.99","size"->entries,"dataset"->++dataset,"type"->"All","analysis"->caption[]}
]
]
];

Monte Carlo: 1 to 4 digits, 0.01-99.99

All the methods above systematically generate all possible numbers in the given range and format, then try all possible ways of keying them. This is slow! An alternative method is to use a Monte Carlo approach. Rather than generate all numbers, we generate numbers randomly, and we try keying them, but making random keying errors with the probabilities p and q. The longer the program is allowed to run, the more accurate its results. To generate numbers, there has to be a known e; it cannot be done symbolically. So we generate data for particular values of e, then use Mathematica to generate interpolation functions in terms of e; this gets us back to the same format of data we have for the purely symbolic approaches — except, instead of an expression in e, we have a expression like F[e], where F is some Mathematica function that fits the data. For either form, we can plot graphs etc easily.

keyError[k_]:=RandomChoice[DeleteCases[{0,1,2,3,4,5,6,7,8,9,dot},k]];

AppendTo[data,
Module[{g,b,di,dj,steps=100,eval=0,prob,i,val,keyPresses,monte,t,experiments=If[fastRun,100,5000]},

Print@ProgressIndicator[Dynamic[monte/(experiments*steps)+eval/steps]];

Print["Generate ",experiments," random numbers, repeating with e = 0 to 1 in steps of ",N[1/steps]];

g=makeArray[Table[0,{steps+1}]];
b=makeArray[Table[0,{steps+1}]];

For[eval=0,eval≤steps,eval++,
Block[{e=eval/steps},
For[monte=1,monte<experiments,monte++,
(* generate random number 0.01-99.99 with 0-2 decimal digits *)
di=PadLeft[IntegerDigits[RandomInteger[{1,9999}]],4];

(* convert to XX.XX format and remove trailing zeros after decimal point; also remove decimal point if two trailing zeros, XX.00 *)
di=If[di[[4]]≠0,{di[[1]],di[[2]],dot,di[[3]],di[[4]]},
If[di[[3]]≠0,{di[[1]],di[[2]],dot,di[[3]]},{di[[1]],di[[2]]},]
];
If[di[[1]]==0,di=Rest[di]];(* delete any leading zero, like 0X.XX *)

(* di is now a list of what user should key; try keying it to make dj *)
keyPresses=Length[di];
val=value[di];
(* now try keying di *)
dj=If[e==0,di,
If[e<1,
Table[RandomChoice[{p,q,ok}->{keyError[di[[i]]],dollar,di[[i]]}],{i,keyPresses}],
Table[RandomChoice[{p,q}->{keyError[di[[i]]],dollar}],{i,keyPresses}],
]];
t=value[dj];(* t is the device's idea of the value of the keyed number *)

(* set t to the scaled out by t ratio so it can be used directly as an index into the bins; if the scaled t is out of range, make it huge *)
t=If[t==0||(t=scale*Max[val/t,t/val])≥huge,huge,Round[t]];

(* if the keyed number is invalid, add to bad bin, else add to good bin *)
If[error[dj],b[[t,eval+1]]+=1,g[[t,eval+1]]+=1]
]
]
];
{"good"->g,"bad"->b,"description"->"Monte Carlo experiments, generating randomly upto 4 digits with or without a decimal point, 0.01-99.99","size"->experiments,"dataset"->++dataset,"type"->"All","analysis"->caption[]}
]
];

Import/export data — generally a lot faster than regenerating it next time!

Save any data that has been worked out. On a first run (or when fastRun is True) this data probably isn't worth saving.

If[!fastRun,
Module[{file=Close[OpenWrite["data.m"]]},
Print["Data is stored in: ",Directory[],"/",file];
Save[file,data];
Print["Data last saved: ",DateString@Date[]]
],
Print["As this is a first run, we didn't bother to save the data! (It might have overwritten data you wanted to keep.)"]
];

Recover data from data file:

If[!fastRun,
Clear[data];
Get["data.m"];
];

See what you have got:

"You have data defined for:\n"<>StringJoin@Map["  * "<>("description"/.#)<>"\n"&,data]

Combine all available results

Since it's tedious to wait for all those results, the code below will combine whichever results have got values, and then average them weighted by how many values each range of data covers.

Data is a vector with elements in this sort of format: {"good" -> goodBins, "bad" -> badBins, "description" -> "e.g., 3 digits with a decimal point after second digit, 10.1-99.9", "size" -> entries, "type" -> "Partial", "dataset" ->5, "analysis"->assumptions}. The goodBins and badBins will contain symbolic expressions in e, or if generated by the Monte Carlo process, they will be vectors of data. If vectors, the process below notices the Head of the data is List and then uses Mathematica's Interpolation function to convert the vectors to functions in e.

At the end of the following code, the cumulative probabilities (i.e., normalized frequencies) are in the vectors prevented and notPrevented and now represented as pairs {x, y}, where y is the probability of an out by x error; note that y is a function of e.

Module[{dg,db,i,j,k,validData=0,steps,types={},otherStuff=""},
Scan[AppendTo[types,"type"/.#]&,data];
If[types=={},
MessageDialog["No data has been generated yet!"];
Return[]
];
If[MemberQ[types,"All"],
(* choose one of them, if any *)
i={"Partial datasets, which will all be combined"->{}}; 
Scan[If["All"==("type"/.#),AppendTo[i,("Dataset "<>ToString["dataset"/.#]<>"\n"<>("analysis"/.#)<>". "<>("description"/.#))->#]]&,data];
i=ChoiceDialog["You have complete "<>If[MemberQ[types,"All"],"and partial ",""]<>"datasets.\n\nWhich do you want to analyze?",i];
If[i≠{},
validData="size"/.i;
{dg,db}={"good"/.i,"bad"/.i}/validData;
];
];
If[validData==0&&MemberQ[types,"Partial"],
Print["You have ",Count[types,"Partial"]," partial data sets, which will be combined."];
(* combine what defined data there are *)
Scan[
If["Partial"==("type"/.#),
Print["Using ","description"/.#," values covering ","size"/.#," samples.\n","analysis"/.#,"."];
t={"good"/.#,"bad"/.#}/("size"/.#);
If[validData>0,{dg,db}+=t,{dg,db}=t];
validData++;
]&,
data];
{dg,db}/=validData
];

If[validData==0,
MessageDialog["You haven't selected any data to analyze!"],
Print[validData," values combined..."];
prevented=makeArray[];
notPrevented=makeArray[];
j=k=0.0;
For[i=huge,i≥1,i--,
j+=db[[i]];
k+=db[[i]]+dg[[i]];
{prevented[[i]],notPrevented[[i]]}=Map[{N[i/scale],#}&,
(* if the Monte Carlo method was used, j and k will be vectors of raw data *)
If[Head[j]===List,
(* values are vectors from Monte Carlo and must be interpolated *)
steps=Length[j]-1;
Interpolation[##,InterpolationOrder->2][e]&/@Map[MapIndexed[Function[{x,y},{N[(y[[1]]-1)/steps],x}],#]&,{j,k}],{j,k}]];
]
];

(* the first few entries correspond to out by factors less than 1, and therefore will be zero, so they are deleted *)
prevented=Drop[prevented,scale-1];
notPrevented=Drop[notPrevented,scale-1];

(* simplification makes everything much faster when we use the formulae later... *)
prevented=Simplify@prevented;
notPrevented=Simplify@notPrevented;
];

(* change the 10 below (both lines) to get data for other sorts of out by r errors *)
good=yintercept[10,prevented][[2]];
bad=yintercept[10,notPrevented][[2]];

This is what the data looks like (see below); you can see it is symbolic, using expressions or functions in e (if it is generated from Monte Carlo data, you will just see the interpolation functions and not the details). I've deleted some tiny constants to make the formulae simpler to read.

Grid[{{Defer[good],"=",approx[good]},{},{Defer[bad],"=",approx[bad]}}]/.approx[number_]:>N[number/.n_?(Abs[#]<(10^(-3))&)->0]//TraditionalForm

We now want to study out by 10 errors in particular, which is what the yintercept[] function gave us a few lines above. If you want to explore out by 5 errors, say, change its parameters above.

Visualize out by 10 results

Red lines in graphs below represent conventional devices; green lines are error-blocking devices (and therefore have better out by 10 performance).

How do we want to visualize results? The ImageSize->1000 parameter in Plot (if not commented out) is useful if you want to cut-and-paste the graphics to (say) PowerPoint, as it makes it bigger and of a consistent size. For presentations, we probably want to provide axes labels in the presentation tool rather than in Mathematica. Hence there is a flag Presentation. Set it to False if you want to live in Mathematica and get axes labels; set it to True if you want to cut-and-paste plots to other applications.

Presentation=ChoiceDialog["Do you want to have labeled graphs, or to have plain graphs at a large, standard size for presentations?",{"Axes labeled, standard size"->False,"No axes labels, large size"->True}];
"Graphs formatted "<>If[Presentation,"large, no axes labels","with labelled axes"]<>"."

Note that you may need to delete plots if you want to resize them — Mathematica usually generates new plots filling the cells the same size as old plots, so if you delete old plots, the new plots will be the size you want them.

Plot[
{good,bad},{e,0,1},
AxesOrigin->{0,0},
PlotStyle->{Directive[Darker[Green],Dashing[{.01,.01}]],Red},
PlotRange->{{0,1},Automatic},
BaseStyle->{Thickness[.005],14},
Evaluate@Sequence[If[Presentation,ImageSize->1000,{PlotLabel->caption[],AxesLabel->{"Keystroke error\nprobability, e","Out by 10\nprobability"}}]]
]

Do the same plot, but now in black-and-white for papers and stuff.

Plot[
{good,bad},{e,0,1},
AxesOrigin->{0,0},
PlotStyle->{Directive[Black,Dashing[{.01,.01}]],Black},
PlotRange->{{0,1},Automatic},
BaseStyle->{Thickness[.005],14},
Evaluate@Sequence[If[Presentation,ImageSize->1000,{PlotLabel->caption[],AxesLabel->{"Keystroke error\nprobability, e","Out by 10\nprobability"}}]]
]

Note. Mathematica does not plot the Monte Carlo model well, as the randomness doesn't handle e→0 very well.

What effective improvement does it give? The graph below shows the saving as e tends to zero; it doesn't work well with Monte Carlo data because of randomness.

Module[{r=good/bad/.e->0.001,xpos=.3},
Plot[good/bad,{e,0,1},
AxesOrigin->{0,0},
PlotStyle->{Darker[Green],Red},
PlotRange->{{0,1},Automatic},
BaseStyle->{Thickness[.005],14},
Epilog->{Thickness[.001],Arrow[{{xpos,r-.2},{0,r}}],Text[" "<>ToString[Round[100r,1]]<>"%",{xpos,r-.2},{-1,0}]},
Evaluate@Sequence[If[Presentation,ImageSize->1000,{PlotLabel->caption[],AxesLabel->{"Keystroke error\nprobability, e","Drop in out by 10 errors\nfor given e"}}
]
]
]
]

For a given value of b (say 0.01) we want to find the e such that bad[e]=good[b], so we can draw a nice graph. Although Mathematica has a Solve function, it doesn't work well on InterpolationFunctions, which result from the Monte Carlo analysis. Here, we use a simple binary search — it's just run often enough to be accurate enough.

BinarySearch[f_,val_]:=
Module[{lo=0,hi=1,h,n},
Do[
h=((f/.e->lo)+(f/.e->hi))/2;
If[h<val,lo=(lo+hi)/2,hi=(lo+hi)/2],
{10}];
N@lo
]

Module[{b=0.1,bb},
bb=BinarySearch[bad,good/.e->b];
Plot[{good,bad},{e,0,2b},
AxesOrigin->{0,0},
PlotStyle->{Directive[Darker[Green],Dashing[{.01,.01}]],Red},
PlotRange->{{0,2b},Automatic},
BaseStyle->{Thickness[.005],14},
Epilog->{
Thickness[.002],
Line[{{b,0},{b,bad/.e->b}}],
Arrow[{{b,bad/.e->b},{b,good/.e->b}}],
Arrow[{{b,good/.e->b},{bb,bad/.e->bb}}],
Line[{{bb,0},{bb,bad/.e->bb}}]
},
Evaluate@Sequence[If[Presentation,ImageSize->1000,{PlotLabel->caption[],AxesLabel->{"Keystroke error\nprobability, e","Out by 10\nprobability"}}
]
]
]
]

Visualize possible improvements

saving=1-good/bad;

monteCarlo[]:=MemberQ[saving,InterpolatingFunction,Infinity,Heads->True]

For the Monte Carlo data, this will give a bumpy plot! (The next plot is a smoothed version.)

Plot[saving,{e,(10^(-3)),1},
AxesOrigin->{0,0},Evaluate@Sequence[If[Presentation,ImageSize->1000,{PlotLabel->caption[],AxesLabel->{"Keystroke error\nprobability, e","Improvement\nblocking errors"}}
]
]
]

The following code smooths the plot, which is very useful for the irregular Monte Carlo data, but probably won't make any difference for symbolic data (unless you didn't generate very much of it).

ListPlot[MovingAverage[Table[saving,{e,(10^(-3)),1,(10^(-2))}],10],Joined->True,AxesOrigin->{0,0},DataRange->{0,1},Evaluate@Sequence[If[Presentation,ImageSize->1000,{PlotLabel->caption[],AxesLabel->{"Keystroke error\nprobability, e","Improvement\nblocking errors"}}
]
]
]

typicalSaving=If[monteCarlo[],
MovingAverage[Table[saving,{e,(10^(-3)),1,0.01}],10][[1]],
saving/.e->(10^(-3))];
CellPrint[If[Abs[typicalSaving]<(10^(-3)),
Cell["Oops, you haven't generated enough data to get good results!", "Subtitle"],Cell["The significance of blocking "<>ToString[Round[100typicalSaving,.1]]<>"% of out by 10 errors...", "Subtitle"]
]
]

Original figures below are taken from Vicente et al (2003).

We convert to per year (divide by 12), convert to UK population (60,943,912), based on US population (303,824,640). The Abbott pump represents (Vicente, correspondence) 75% of the US market, so we could legitimately scale by 4/3 to get representative rate for all PCA devices. (NB I don't believe these population counts to ±1 person, but they'll do the job!)

If typing error rates are causing a rate of saving bad doses that can be detected,  then the detection rate will eliminate some fatalities:

italicize[n_List]:=Map[italicize,n];
italicize[n_]:=If[StringQ[n],Style[n,Italic],n];

embolden[n_List]:=Map[embolden,n];
embolden[n_]:=If[StringQ[n],Style[n,Bold,20,FontFamily->"Helvetica"],n];

display[d_,title_]:=
Print@Text@Grid[{
{embolden@title,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft},
italicize@{"","","Estimates of true incidence",,"Estimates of prob",SpanFromLeft},
italicize@{"","Reported deaths","Low","High","Low","High"},
Prepend[d[[1]],italicize@"Min"],
Prepend[d[[2]],italicize@"Max"]},
Alignment->{{Center,Right,Right,Right,Right,Right},Baseline},Dividers->{All,{1->True,3->True,-1->True}}];

vicente=
{
{5,65,417,2.95 10^-6,1.89 10^-5},
{8,104,667,4.72 10^-16,3.03 10^-5}
};

display[vicente,"Vicente original data, covering 12 years in USA"];
display[vicentePA=Map[Round[#*{c,c,c,1,1}/.c->(1/12)/303.824640,.001]&,vicente],"Fatalities, scaled per 1M population pa for one PCA"];
display[Map[Round[#*{typicalSaving,typicalSaving,typicalSaving,0,0},.0001]&,vicentePA],"Prevented fatalities per 1M pa for one PCA device\n"<>caption[]];display[Map[Round[#*{typicalSaving,typicalSaving,typicalSaving,0,0},.0001]&,60.943912*vicentePA],"Prevented fatalities UK pa for one PCA device\n"<>caption[]];

Visualizing the data another way ...

zapNullOptions[] takes a list of option values, and is a handy way to clear empty options.

zapNullOptions[n_]:=DeleteCases[n,False->_]/.(True->arg_)->arg;

explainPlot[], below, draws a sequence of graphs to help explain the data. Useful for PowerPoint talks!

explainPlot[eValue_,detail_]:=
Module[{ratio=10,icpt,nppt,fig,explain,saving,preventedate,notpreventedate},
preventedate=prevented/.e->eValue;
notpreventedate=notPrevented/.e->eValue;
fig[n_]:=ToString[Round[n,0.001]];
explain[t_,xy_]:=
Module[{dh=1},
{Text[t,xy+{-0.1,-0.00075},{1,1}],Arrow[{xy,{0,xy[[2]]}}]}
];
icpt=yintercept[ratio,preventedate];
nppt=yintercept[ratio,notpreventedate];
saving=(nppt[[2]]-icpt[[2]])/nppt[[2]];
Print@ListPlot[zapNullOptions@
{{{1,0},{1,1}},
detail>1->{{ratio,0},{ratio,1}},
detail>2->preventedate,
notpreventedate
},
Evaluate@Sequence[If[Presentation,ImageSize->1000,{PlotLabel->caption[],AxesLabel->{"Out by x","Probability of out by x"}}]],
Joined->True,
PlotStyle->zapNullOptions@
{{Thickness[.0025],Dashing[{.01,.01}],Gray},
detail>1->{Thickness[.0025],Dashing[{.01,.01}],Gray},
detail>2->Darker[Green],
Red},
PlotRange->{{0,11},{0,0.022}},
BaseStyle->{Thickness[.005],14},
AxesOrigin->{0.00,0},
Ticks->{Range[10],Automatic},
AspectRatio->.5,
Epilog->zapNullOptions@
{Thickness[.002],Arrowheads[.02],
detail>3->explain["Safer device, Prob(out by "<>fig[ratio]<>")="<>fig[icpt[[2]]],icpt],
detail>1->explain["Ordinary device, Prob(out by "<>fig[ratio]<>")="<>fig[nppt[[2]]],nppt]
}];
If[detail==4,CellPrint[Cell[
"For e="<>fig[eValue]<>", with prob "<>fig[nppt[[2]]]<>" a normal device is out by "<>ToString[ratio]<>"; with improved user interface with error prevention it is probability "<>fig[icpt[[2]]]<>" it is out by "<>fig[ratio]<>" That is, for a ratio "<>fig[ratio]<>", it is "<>fig[100saving]<>"% safer.", "Text"]]]
];

If[!fastRun,explainPlot[0.01,#]&/@Range[4];,"You wanted quick results; these graphs aren't very informative on small datasets, so have not been drawn."]

References

K. J. Vicente, K. Kada-Bekhaled, G. Hillel, A. Cassano & B. A. Orser, "Programming errors contribute to death from patient controlled analgesia: Case report and estimate of probability," Canadian Journal of Anesthesiology, 50(4):328–332, 2003.

K. J.  Vicente, Correspondence, Canadian Journal of Anesthesiology, 50(8):856-857, 2003.