\newpage
\appendix
This appendix describes our ``STAR2'' AI tool that we use to generate conclusions
from these models.
A simulated annealer is used to  sample that space and a back select
algorithm is used to  rank feature ranges according to how well they reduce effort/ defects/ 
threats. This section will offer a brief explanation of these STAR2 internals.

{\em Monte Carlo} algorithms randomly sample the space
of possible controllable model states. 
A  {\em Metropolis} Monte Carlo algorithm~\cite{metropolis53}
creates new states
by  small mutations to some
{\em current} state. 
If a new state is ``better'' (as assessed via an
{\em energy function}),  it becomes the new {\em current}
state used  for future
mutations.  Otherwise, a Boltzmann acceptance
criteria is used to probabilistically   decide to
assess the new state: the worse the new state, the less
likely that it becomes the new current state.
The algorithm is silent on the mutation mechanism. For our experiments,
we freeze $\frac{2}{3}$ of the features and
randomly  select ranges for the rest.

In 1983, Kirkpatrick et.el.~\cite{kirkpatrick83} 
proposed a modification that was
inspired
by a procedure used to make the strongest possible glass. Initially,
glass is heated to allow atoms to move freely. The temperature
is then slowly lowered such that the atoms can find the stablest
configuration with lowest energy. 
A {\em simulated annealing} (SA) algorithm adds
a ``temperature'' variable to the Boltzmann accept criteria
such that, at high temperatures, it is more likely that the
algorithm will jump to a  new worst current state.
This allows the algorithm to jump out of local minima while
sampling the space of options. As the temperature cools,
such jumps become less likely and the algorithm reverts
to a simple hill climber.

With SA, better solutions have lower ``energy''
so the fitness function $E$ is defined such that $E \ge 0$ and
{\em lower} values are
$better$. 
For our purposes, we model effort ($Ef$) as the distance to
the origin of the 4-D space with defects ($De$), threats ($Th$), and schedule ($Mo$):
\begin{equation}\small\label{eq:eng}
\begin{array}{c}
RD = relydefect^{(RELY-3)}\\\\
E=\frac{\left(\sqrt{(\overline{Ef}*\alpha)^2 + (\overline{De}*(\beta+RD))^2 + (\overline{Th}*\gamma)^2 + (\overline{Mo}*\delta)^2}\right)}{\sqrt{\alpha^2 + (\beta+RD)^2 + \gamma^2 + \delta^2}}\\
\end{array}
\end{equation}
Here, $\overline{x}$ is a normalized value 
$0 \le \frac{x - min(x)}{max(x) - min(x)} \le 1$. 
Hence, our energy ranges $0 \le E \le 1$ and
{\em lower} energies are {\em better}. $\alpha$, $\beta$, $\gamma$, and $\delta$ 
are the weighting factors for the different models. $relydefect$ is a constant
that defines a relationship between RELY and the defect model: as RELY increases,
the weighting of the defect model increases as well.

The acceptance criteria $P$ for a new state is defined using
the current temperature $T$, 
the energy of the current solution ($e$),
and the energy of the new mutation ($en$):
\begin{equation}\small
P(e,en,T) = e^{(e -en)/T}
\end{equation}
$T$ is defined to $decrease$ as the simulator loops through $k= 1\ldots kmax$ iterations.
We use
$T = e^{-100*k/kmax}$.

\begin{figure}[!t]
\scriptsize
\begin{verbatim}
function sa(kmax) 
  s := s0; e := E(s)           // Initial state, energy.
  sb := s; eb := e             // Initial "best" solution
  k := 0                       // Energy evaluation count.
  while k < kmax and e > emax  // Loop
    sn := neighbour(s)         //   Pick some neighbour.
    en := E(sn)                //   Compute its energy.
    if en < eb then            //   Is this a new best?
      sb := sn; eb := en       //     Yes, save it.
    if random() < P(e, en, temp(k/kmax)) 
     then s := sn; e := en     //     Maybe jump
    k := k + 1                 //   One more evaluation done
  return sb                    // Return best
\end{verbatim}
\caption[SA pseudo-code]{SA pseudo-code:
a new solution $sn$ (with new energy $en$) replaces the current solution                                  
if (a)~it has a lower energy or (b)~the acceptance predicate $P$
endorses it. Only in the case of (a)
should the new solution
replaces the current best solution.
} \label{fig:sa}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=2.5in]{jpl1.pdf}
%eps/JPL.pdf}
\end{center}
\caption{Processing a JPL requirements model. }\label{fig:jpl}
\end{figure}
Two advantages of SA are their {\em implementation simplicity} and their {\em ability to handle
non-linear models}:
\bi
\item {\em Implementation simplicity:}
\fig{sa} illustrates the simplicity of the algorithm. Memory is only required
for one current solution ($s$), one new solution ($sn$) and one best solution ($sb$) that stores
the best solution seen at any time in the simulation.
\item
{\em Non-linear models:}
 Previously~\cite{fea02a,me03i}, we have applied SA to non-linear JPL requirements models
where minimizing the cost of project mitigations can decrease the number of requirements
achieved by that project. Hence, decreasing {\em both} the cost and achieved requirements
is a non-linear problem that must trade between minimizing cost and increasing
requirements coverage.
The top-left line  of \fig{jpl} divides the behavior of the JPL requirements models
before  and after  simulated annealing. 
As shown below the line, 
initial Monte Carlo sampling of the possible mitigations lead to a large range
of costs and benefits. Simulated annealing found a set of mitigations
that lead to the small cloud of solutions above the line. Compared to the initial samples, these
new solutions had decreased cost, increased benefit (number of requirements covered), and decreased variance (shrank
the space of solutions). 
\ei
Two disadvantages of SA algorithms are their {\em incompleteness} and the {\em
complexity}
of their solutions:
\bi
\item 
{\em Incompleteness:} In our domain, we have some evidence that the incomplete nature of the heuristic SA
search is not a major problem.
\fig{saeg} shows a sample run of our SA tool running on our prediction models
for $K=10,000$ simulations. 
As $k$ increases for $1{\ldots}10,000$, it becomes less and less
likely that a better $best$ has been missed. Hence, we run our simulations
for ten times the period it takes for $best$ to stabilize (at $k\approx 1000$). 
\begin{figure}
\begin{center}
\includegraphics[width=2.7in]{saeg1.pdf}
\end{center}
\caption[SA solutions]{Dots \& lines are
SA output from current \& best solution (respectively)  after $k$ simulations.
}\label{fig:saeg}
\end{figure}
\item
{\em Solution complexity:}
Simulated annealers offer constraints to all controllable features. 
Often this
is an over-constrained solution since, 
in many domains, a repeated empirical result is a {\em feature subset selection}
effect; i.e. models that constrain M variables perform just as well, or better,
than models that constrain N variables ($M\ll N$)~\cite{kohavi97,hall03}.
For example, Kohavi~\cite{kohavi97}
studied some machine learners to find  that using just 19\% of the available features
{\em increased} prediction accuracy by just 2.14\% (on average).
For another example,
when feature subset selection was applied to the JPL requirements model of \fig{jpl}, we found that
up to $\frac{2}{3}$-rds of the features can be left unconstrained,
without effecting the conclusions~\cite{fea02a}. 
\ei
%
%\begin{figure}[!t]
%\scriptsize
%\begin{center}
%\begin{tabular}{|l|r|r|r|p{0.95cm}|}\cline{2-4}
%\multicolumn{1}{c}{}     & \multicolumn{3}{|c|}{av.  number of features} \\\hline    
%dataset         &before &after &$\frac{after}{before}\%$ &accuracy change\\\hline    
%breast cancer &10  &2.9 &29\% &+0.14\%   \\\hline    
%cleve           &13  &2.6 &2\%  &+5.89\%   \\\hline 
%   crx             &15  &2.9 &19\% &+4.49\%   \\\hline  
%  DNA             &180 &11  &6\% &+3.63\% \\\hline    
%horse-colic     &22  &2.8 &13\% &+1.63\%   \\\hline
% Pima            &8   &1   &13\% &+0.79\%   \\\hline   
% sick-euthyroid &25  &4   &16\% &+0.38\%   \\\hline    
%soybean         &35  &12.7 &36\% &+0.15\%  \\\hline\hline    
%average         &38.5&4.99 &19\% &+2.14\% \\\hline
%\end{tabular}
%\end{center}
%\caption{Feature subset selection results from~\cite{kohavi97}}\label{fig:fssr}
%\end{figure}

In terms of  the goal of this paper, exploring feature subset selection
is very important.  Before we can {\em remove} non-essential features,
we must first  {\em rank} them according to their effectiveness.
As we shall see, this ranking will be insightful to the task of
assessing the relative value of new
 ASE
tools and other methods
for reducing effort/ defects/ threats.


\textbf{There are two ranking methods that can be used by STAR2. These are {\em Bayesian ranking} and {\em Energy ranking}. The former is the main ranking method used by STAR2, while the latter is included in STAR2 for verifying that the Bayesian ranking is valid.}


STAR2 ranks the features ranges
seen in $K$ runs of a simulated annealer by dividing the $K$ runs into:
\bi
\item {\em Best:} those associated with the  BEST\% solutions
(i.e. those with the BEST\% {\em least} energy); 
\item And the {\em rest} (i.e. the other 100-BEST\% of solutions).
\ei
It then computes the probability that a range is found in $best$.
using Bayes' Theorem.
Informally,
the
theorem says $next=old*new$ i.e. what we'll believe {\em next}
comes
from how {\em new} evidence effects {\em old} beliefs.  More
formally:\begin{equation}\label{eq:n}\small
P(H|E)=P(E|H) P(H)\; /\; P(E)\end{equation}
i.e. using
evidence $E$ and a prior probability $P(H)$ for
hypothesis $H\in\{best, rest\}$.
The theorem calculate a  posteriori  probability
$P(H|E)$. 
Simple  Bayes classifiers are often called ``n\aive'' since they assume
independence of each feature. While this assumption
 simplifies the implementation  
(frequency
counts are required only for each feature),
 it is possible
that
correlated events are missed by this ``n\aive'' approach.
Domingos and Pazzani show theoretically
that the independence assumption is a problem in a vanishingly small percent
of cases~\cite{domingos97optimality}. This explains the repeated
empirical result that, on average, seemingly n{\aive} Bayes classifiers perform as well
as other seemingly more sophisticated schemes (e.g. see Table~1 in~\cite{domingos97optimality}).


When applying the theorem, {\em likelihoods}
are computed from observed
frequencies, then normalize to create
probabilities (this normalization cancels out $P(E)$ in \eq{n}, so it need not be computed).
For example, after $K=10,000$ runs divided into  1,000 lowest 10\% {\em best} solutions and 9,000
{\em rest}, the  range $rely=vh$ might appears 10 times in the $best$ solutions, but only
5 times in the $rest$. Hence:

{\scriptsize
\begin{eqnarray}\nonumber
E&=& (reply=vh)\\\nonumber
P(best)&=&1000/10000 = 0.1\\\nonumber
P(rest)&=&9000/10000 = 0.9\\\nonumber
freq(E|best)&=&10/1000 = 0.01\\\nonumber
freq(E|rest)&=&5/9000 = 0.00056\\\nonumber
like(best|E)&=&freq(E|best) \cdot P(best)=0.001\\\nonumber
like(rest|E)&=&freq(E|rest) \cdot P(rest)= 0.000504\\\label{eq:ten}
P(best|E)&=& \frac{like(best|E)}{like(best|E) + like(rest|E)}=0.66
\end{eqnarray}}
Previously~\cite{clark05}
we have found that  \eq{ten} 
is a poor ranking heuristic
since it is distracted by low frequency evidence.
For example, note how
the probability of $E$ belonging to the best class
is moderately high even though its support is very low; i.e. $P(best|E)=0.66$ but
$freq(E|best) = 0.01$. 

To avoid such unreliable low frequency evidence, we augment \eq{ten}
with a support term.
Support should $increase$ as the frequency of a range $increases$,  i.e. 
$like(x|best)$ is a valid support measure.
STAR2 hence ranks ranges via
\begin{equation}\label{eq:ps}\small
P(best|E) * support(best|E) = \frac{like(x|best)^2}{like(x|best) + like(x|rest)}
\end{equation}


While the main ranking algorithm used by STAR2 is the Bayesian method 
explained above, another very basic and intuitive baseline method was 
implemented to test Bayesian ranking. This method is Score-Based Energy
Ranking. Basically this ranking scheme uses the means of the results 
from the Simulated Annealing simulations to assess and rank COCOMO II 
attribute-value combinations. The mean results are based the scores 
generated by the Energy Function of the SA algorithm. These scores are 
basically the distance from the origin to a point in three dimensional 
space whose coordinates are the normalized Effort, Defects, and Threat scores from our models. Energy Ranking works as follows:

\bi
\item[1.] {\em DATA COLLECTION:}
This is done while the SA algorithm is 
running whenever the algorithm generates a new candidate best result.
 The attribute-value combinations that are used for that candidate are 
detected, and data concerning the score is stored for each of these used attribute-value combinations.
\newline
\item[2.] {\em MEAN SCORE CALCULATION:}
After the SA goes through its K 
runs, the mean scores for all the attribute-value combinations are calculated. 
After this step, the attribute-value combinations are ready to be ranked.
\ei


While Energy Ranking seems to be simple enough, is has one main 
potential weakness: outliers in simulations could taint the score of
an otherwise good attribute-value combination, causing its rank to 
decrease compared to seemingly worse combinations. This issue is not
present in Bayesian Ranking since it is dependent on frequency and support rather than scores.

To apply \eq{ps}, STAR2 runs in six phases.
In terms of standard machine 
learning theory, step 1 generates a training set; steps 2,3,4 do some generalization;
and step 5 tests the learned theory on data not seen during training.

\begin{description}
\item[1.] {\em SAMPLE:}
To sample the ranges from the models,
STAR2 runs the simulated annealer  $K_1$ times.
\item[2.] {\em DISCRETIZE:}
The data seen in the $K_1$ samples is then discretized into $D=10$ bins.
Discretization converts a continuous range into a histogram
with $n$ break points $b_1\ldots b_n$ where $\left(\forall i < j : b_i \le b_j\right)$. 
After discretization,
many observations can fall into the same range
between $b_i$ and $b_{i+1}$
at frequency counts $c_i$. This study used equal width discretization; i.e. 
\[\forall i,j:
\left(b_i - b_{i-1}\right) = 
\left(b_j - b_{j-1}\right)\]
\item[3.] {\em CLASSIFY:}
The ranges are then classified
into those seen in BEST\% {\em best} or {\em rest}.
\item[4.]
{\em RANK:}
The ranges are then ranked in increasing order either using Support-Based Bayesian Ranking \tion{sec:bayrank}  or Score-Based Energy Ranking \tion{sec:enrank}.
\item[5.]
{\em PRUNE:}
STAR2 then runs $K_2$ experiments with the models where the top ranked ranges
$1..X$ ranges are pre-set and the
remaining ranges can be selected at random.
In terms of the business case for automated software engineering tools,
this step is crucial since it would rank new ASE tools alongside 
other methods. 
\item[6.]
{\em REPORT:}
STAR2 returns the $1..X$ ranges that  produce the least effort, defects, and threats.
\end{description}

To run our experiments, we had to apply our engineering judgment to set the parameters. The following are the default values:
\[K_1=10,000, K_2=1,000, D=10, BEST=10\%\]

%Initially, we planned experiments that varied these parameters,
%as well as trying other discretization policies.
%However, our initial results were so promising (see below) that we are not currently motivated
%to explore 
%other parameter settings (such an exploration is planned for future work).

Note that for Energy Ranking, the algorithm is slightly different. Energy ranking does not require the disretizing and the classifying steps 


%\begin{figure}[!t]
%\begin{center}
%\scriptsize
%\begin{tabular}{c|lrr|lr}
%       &\multicolumn{3}{c|}{ranges}      &\multicolumn{2}{c}{values}\\\hline
%project&feature&low&high&feature&setting\\\hline
%&rely&3&5&tool&2\\
%Flight:&data&2&3&sced&3\\
%&cplx&3&6&&\\
%&time&3&4&&\\
%&stor&3&4&&\\
%&pvol&2&4&&\\
%&acap&3&5&&\\
%&apex&2&5&&\\
%&pcap&3&5&&\\
%&plex&1&4&&\\
%&ltex&1&4&&\\
%&pmat&2&3&&\\
%&Ksloc&7&418\\\hline
%&rely&1&4&tool&2\\
%Ground:&data&2&3&sced&3\\
%&cplx&1&4&&\\
%&time&3&4&&\\
%&stor&3&4&&\\
%&pvol&2&4&&\\
%&acap&3&5&&\\
%&apex&3&5&&\\
%&pcap&3&5&&\\
%&plex&1&4&&\\
%&ltex&1&4&&\\
%&pmat&2&3&&\\
%&Ksloc&11&392\\\hline
%&prec&1&2&data&3\\
%OSP:&flex&2&5&pvol&2\\
%Orbital&resl&1&3&rely&5\\
%space&team&2&3&pcap&3\\
%plane&pmat&1&4&plex&3\\
%&stor&3&5&site&3\\
%&ruse&2&4&&\\
%&docu&2&4&&\\
%&acap&2&3&&\\
%&pcon&2&3&&\\
%&apex&2&3&&\\
%&ltex&2&4&&\\
%&tool&2&3&&\\
%&sced&1&3&&\\
%&cplx&5&6\\
%&KSLOC&75&125\\\hline
%&prec&3&5&flex&3\\
%OSP2 &pmat&4&5&resl&4\\
%&docu&3&4&team&3\\
%&ltex&2&5&time&3\\
%&sced&2&4&stor&3\\
%&KSLOC&75&125&data&4\\
%&&&&pvol&3\\
%&&&&ruse&4\\
%&&&&rely&5\\
%&&&&acap&4\\
%&&&&pcap&3\\
%&&&&pcon&3\\
%&&&&apex&4\\
%&&&&plex&4\\
%&&&&tool&5\\
%&&&&cplx&4\\
%&&&&site&6\\
%\end{tabular}
%\end{center}
%\caption{Four of the five case studies.}\label{fig:cases}
%\end{figure}



