# Automatic Assessment of the Design Quality of Python Programs with
Personalized Feedback

J. Walker Orr

George Fox University<EMAIL_ADDRESS>Nathaniel Russell

George Fox University<EMAIL_ADDRESS>
###### Abstract

The assessment of program functionality can generally be accomplished with
straight-forward unit tests. However, assessing the design quality of a
program is a much more difficult and nuanced problem. Design quality is an
important consideration since it affects the readability and maintainability
of programs. Assessing design quality and giving personalized feedback is very
time consuming task for instructors and teaching assistants. This limits the
scale of giving personalized feedback to small class settings. Further, design
quality is nuanced and is difficult to concisely express as a set of rules.
For these reasons, we propose a neural network model to both automatically
assess the design of a program and provide personalized feedback to guide
students on how to make corrections. The model’s effectiveness is evaluated on
a corpus of student programs written in Python. The model has an accuracy rate
from 83.67% to 94.27%, depending on the dataset, when predicting design scores
as compared to historical instructor assessment. Finally, we present a study
where students tried to improve the design of their programs based on the
personalized feedback produced by the model. Students who participated in the
study improved their program design scores by 19.58%.

###### keywords:

Assessment, neural networks, intelligent tutoring

##  1 Introduction

Recently there has a been a lot of work in the development of tools for
education in programming and computer science. Specifically there are many
systems for intelligent tutoring which are designed to help students learn how
to solve a programming challenge. The tutoring involved is primarily focused
in suggesting functional improvements, that is, how to finish the program so
that it works correctly.

Intelligent tutors such as [4] uses reinforcement learning to predict a useful
hint in the form of an edit to a student’s program that will get them one step
closer to the goal of a functioning program. It uses histories of edits made
by students, starting with a blank slate and ultimately terminating with a
functional program to train the model. The system is based on Continuous Hint
Factory [9] which uses a regression function to predict a vector that
represents the best hint then translates that vector into a human-readable
edit. Similarly [11] used a neural network to embed programs and predict the
program output. Using that model of the program output, an algorithm was
developed to provide feedback to the student on how to correct their program.
Also, [16] use a recurrent neural network to predict student success at a task
given a history of student submissions of their program for evaluations.

All these systems model student programs from Hour of Code [2]. Hour of Code
is a massively open online course platform that teaches people how to code
with a visual programming language. The language is simple and does not
contain control constructs such as loops.

Moreover, the combination of language and problem setting are simple enough
that there is a single or very few functional solutions for each problem [4].
This level of simplicity precludes the consideration of program design.
However for general purpose programming languages such as Python, there are
many ways of creating functionally equivalent programs. It is important for
the sake of maintainability, modularity, clarity, and re-usability that
students learn how to design programs well.

When it comes to the quality of design, there are varying standards. Further,
some standards are more objective or easier to precisely identify that others.
For example, the use of global variables are both widely recognized as poor
design and are easy to identify. For some programming languages, “linters”
exist to apply rules to check for common design flaws. For Python, Pylint [12]
is a code analysis tool to detect common violations of good software design.
It detects design problems such as the use of global variables, functions that
are too long or take too many arguments, and functions that use too many
variables. Pylint is design to enforce the official standards of the Python
programming community codified in PEP 8 [15].

There are aspects of good design that are difficult to identify. For example,
simple logic is a good design idea, but is quite nebulous. The complexity of a
program’s logic is contextual, it entirely depends on the problem the program
is solving. Also, modularity is universally judged as a quality of good
design, however it is not always clear to what extent a program should be made
modular. How many functions or classes are too many? Again, it depends on the
context of the problem for which the program is designed.

In a professional setting, code reviews are often practiced to promote quality
design that goes beyond the straight-forward rules of “linters.” Code reviews
are a manual process which require a lot of human effort. A recently developed
system call DeepCodeReviewer [5] automates the code review process with a deep
learning model. By using proprietary data on historical code reviews taken
from a Microsoft software version control system, DeepCodeReviewer was trained
to successfully annotate segments of C# code with useful comments on the
code’s quality.

However, to our knowledge, there is no system to perform in-depth code
analysis for the purposes of evaluating and assessing design for general
purpose languages in an educational context. The process of assessing the
design of a program is time consuming for instructors and teaching assistants
and it is an important component of complete intelligent tutoring system. Such
a system needs to be adjusted or calibrated for the context of particular
problems or assignments since there are important aspects of software design
are context dependent. Moreover the system needs to match the particular
standards of an instructor. Hence we propose a system that models design
quality with a neural network trained on previously assessed programs.

###  1.1 Our Approach and Contribution

We propose a design quality assessment system based on a feed-forward neural
network that utilizes an abstract syntax tree (AST) to represent programs. The
neural network is a regression model that is trained on assessed student
programs to predict a score between zero and one. Each feature the model uses
is designed to be meaningful to human interpretation and is based on
statistics collected from the program’s AST. We intentionally do not use deep
learning as it would make the representation of the program difficult to
understand. Personalized feedback is generated based on each feature of an
individual program. By swapping a feature’s value for an individual program
with the average feature value of good programs, it is possible to determine
which changes need to be made to the program to improve is design. The primary
contributions of this work are the following:

  * •

The first to explicitly predict the design quality of programs in an
educational setting to the best of our knowledge.

  * •

High efficacy with an accuracy from 83.67% to 94.27% with only small amounts
of training data required.

  * •

The first intelligent tutoring system for design quality for Python.

  * •

Personalized feedback without the explicit training or annotation.

##  2 Method

Figure 1: The model of program design quality, a feed-forward neural network.
The “Input Layer” is the feature vector created from the AST. The “Hidden
Layer” corresponds to calculation of $x^{\prime}$ specified in Equation 1.
Finally, the “Output Layer” produces a single value, the design score as found
in Equation 2.

The task is to predict a design quality score for a student program written in
Python. The score $y$ is a real number between zero and one. The program is
represented by a feature vector $\vec{x}$ produced by the output of a series
of feature functions computed from the program’s AST.

For an AST $T$, a series of feature functions $f_{i}((T)$ output is
concatenated in to a feature vector $\vec{x}$ that represents key aspects of
the program’s design. The model $g(\vec{x};\Theta)$ is a feed-forward neural
network with a single hidden layer. It is a regression model that predicts the
score $y$ based on the feature vector $\vec{x}$ and parameters $\Theta$.

###  2.1 Features

Figure 2: An abstract syntax tree for a segment of Python code that computes
the average of the values in a list. The AST has been condensed for the sake
of brevity and the restrictions of space. The related code segment is shown in
blue.

Despite recent advances in deep learning, we chose to represent the student
program with feature functions computed on its AST. Deep learning is highly
effective at learning useful feature representations of everything from images
to time series to natural language texts. However, deep learning also requires
large amounts of data and in this setting the quantity of manual annotated
student programs is limited.

Additionally, AST are a natural and effective means of representing and
understanding programs and can be created with free, available tools. An AST
is an exact representation of the source code of program based on the
programming language’s grammatical structure. Producing an AST representation
of a programming language is an essential first step in compilers and
interpreters. The AST of a program contains all the content of its source but
also is augmented with the syntactic relationships between every element. A
parser and tokenizer to produce an AST for the Python programming language is
provided by its own standard library. This makes the AST the natural
representation to use, since it is free, convenient, exact, interpretable, and
does not required any additional data. In contrast, deep learning would
require a large amount of data to effectively reproduce the same
representation.

Prior to representation as an AST, a program must first be broken into a
series of tokens via the process of lexicalization. Lexicalization is the
process of reading a program, character-by-character and dividing into work-
like tokens. These tokens are also assigned a type such as a function call or
variable reference. The grammatical rules of the language are applied to the
lexicalized program to create the AST. In an AST, the leaf nodes of the tree
are the program’s tokens while the interior nodes correspond to syntactic
elements and constructions. For example, an interior node could represent the
body of a function or the assignment of a value to a variable. An example of
an AST is found in Figure 2.

Given that an AST is a complete representation of a program, it is a natural
basis for assessing the quality of a program’s design. Deep learning may be
able to automatically learning the same key syntactic relationships with
enough data, however this information is simply available via AST. Further,
features computed from the AST will be human interpretable unlike a
representation produced by deep learning.

The features we created are all based on statistics collected from a program’s
AST. Some consist of simply counting the number of nodes of a given type, for
example, the number of user defined functions. Other feature functions are
based of subsections of the AST, such as the number of nodes per line or per
function. Finally, some features are ratios or percentages such as the average
percent of lines in a number in a function that are empty. All of the features
are relatively simple and fast to compute, yet generally capture the design
and quality of a program. Each of the feature functions $f_{i}(T)$ we defined
can be found in Appendix A.

###  2.2 Model

The model is a feed-forward neural network [13] with a single hidden layer and
single neuron in the output layer. The model’s structure is illustrated in
Figure 1. The values of the input layer are the feature vector $\vec{x}$. Each
neuron in the hidden layer $x^{\prime}_{j}$ defined with the following
equation:

$\displaystyle
x^{\prime}_{j}=\text{ReLU}\bigg{(}\sum\limits_{i=1}^{d}w_{i,j}x_{i}\bigg{)}$
(1)

where $d$ is the dimension of $\vec{x}$ and $w_{i,j}\in\Theta$ are the
parameters, “weights” of the neuron. We use the ReLU [8] as the activation
function for the hidden layer neurons. The final prediction of the design
score is made by the output layer’s single neuron:

$\displaystyle
y=\sigma\bigg{(}\sum\limits_{j=1}^{d^{\prime}}w_{j}x^{\prime}_{j}\bigg{)}$ (2)

where $d^{\prime}$ is the number of hidden layer neurons and $w_{j}\in\Theta$
are weights of this neuron. The function sigmoid is used because its domain
spans from $(-\infty,\infty)$ but its range is $[0,1]$ which ultimately
guarantees the model always outputs a valid score. The model is trained with
mean squared error as the loss function:

$\displaystyle\text{MSE}=\frac{1}{n}\sum\limits_{k=1}^{n}(y^{*}_{k}-y_{k})^{2}$
(3)

where $y^{*}_{k}$ is the ground-truth design score for the $k^{\text{th}}$
instance i.e. program and $n$ is the number of instances in the training data.
The model is trained with the ADAM algorithm [7] and each parameter in the
model was regularized according to their L2 norm [6]. For all our experiments,
a hidden layer of size 32 was used. The model was trained for 250 epochs and
the model from the best round according to a development set was selected for
our experiments.

###  2.3 Ensemble

Due to the fact that fitting neural networks to data is a local optimization
problem, the effect of initial values of the parameters $\Theta$ of the model
remain after training. The process of training a neural network will produce a
different model given the same data. This variation in the results of a
trained model is particularly pronounced when the training data set is
relatively small. To address this variation and mitigate its impact an
ensemble of models can trained, each with different initial parameter values.
Each model is independently trained and a single prediction is made by a
simple of the average of individual predictions i.e.

$\displaystyle y=\frac{1}{m}\sum\limits_{l=1}^{m}y_{l}$ (4)

where $m$ is the number of models and $y_{l}$ is the prediction of the
$l^{th}$ model. For our experiments, an ensemble of 10 models was used.

###  2.4 Personalized Feedback

The goal of intelligent tutors is to provide personalized feedback and
suggestions on how to improve a program. The most straight-forward means of
providing feedback would be to simply predict which possible improvements
apply to a given program. However, training a model to directly predict
relevant feedback would require a dataset of program with corresponding
feedback and such a dataset can be hard to find or is expense to construct.

In order to avoid the need for a dataset with explicit feedback annotation, we
use our model, trained on predicting design score, to evaluate how changes in
program features would lead to a higher assessed score. Using the training
data, we compute an average feature vector $\underline{\vec{x}}$ of all the
“good” programs i.e. those with a design score greater than 0.75. To generate
feedback for a program, its feature vector $\vec{x}$ is compared to the
average $\underline{\vec{x}}$. For each feature, a new vector
$\vec{x}^{\prime}$ is created by replacing the feature value $x_{i}$ with
value with the average’s value $\underline{x}_{i}$. This process is setting up
a hypothesis, what if the program was closer to the average “good” program
with regards to a particular feature? To answer this, the trained model
$g(\vec{x};\vec{\theta})$ is used to predict a design score for the new vector
i.e. $y^{\prime}_{i}=g(\vec{x}^{\prime};\vec{\theta})$. By comparing the
original score of the program $y$ with the new score $y^{\prime}_{i}$, the
hypothesis can be tested. If the new score $y^{\prime}_{i}$ is greater than
the original predicted score $y$, then the alteration of $x_{i}$ to be closer
to $\underline{x}_{i}$ is an improvement. Feedback based on this alteration is
recommended to the student as personalized feedback. Since each feature in
$\vec{x}$ is understandable to a human, feedback is given in the form of the
suggestion to increase or decrease particular features. The suggestion for
alteration is based on the comparison of $x_{i}$ versus $\underline{x}_{i}$,
if $x_{i}>\underline{x}_{i}$, the feedback of decrease $x_{i}$ is given. In
the other case, where $x_{i}<\underline{x}_{i}$ the feedback is to increase
$x_{i}$. Based on the feature and the feedback of increase or decrease, a
user-friendly sentence is selected from a table of predefined responses. For
example, if $x_{i}$ is the number of user defined functions and $x_{i}=3$,
$\underline{x}=5$, and $y^{\prime}_{i}>y$ then the feedback of “increase the
number of user defined functions” is created.

##  3 Experiments

Method | Travel | Budget | RPS | Craps | Combined
---|---|---|---|---|---
| MSE | Accuracy | MSE | Accuracy | MSE | Accuracy | MSE | Accuracy | MSE | Accuracy
Linear Regression | 0.009 | 93.09% | 0.032 | 87.43% | 0.038 | 83.2% | 0.043 | 84.06% | 0.027 | 87.03%
Decision Tree | 0.018 | 90.10% | 0.031 | 87.26% | 0.078 | 77.33% | 0.072 | 79.41% | 0.076 | 81.42%
Sig. Linear Regression | 0.022 | 89.60% | 0.046 | 85.13% | 0.086 | 77.91% | 0.063 | 81.43% | 0.070 | 80.64%
Neural Network | 0.007 | 93.48% | 0.024 | 88.48% | 0.041 | 83.9% | 0.08 | 79.57% | 0.033 | 85.61%
Ensemble | 0.005 | 94.27% | 0.022 | 90.14% | 0.022 | 87.66% | 0.053 | 83.67% | 0.031 | 86.99%
Table 1: Design Score Prediction Results

The system was evaluated in two different experimental settings. The first
evaluation is direct test of the model’s accuracy on known design scores. For
this, several different datasets and settings were compared against several
baselines. The second evaluation is a small study of how student’s responded
to the system’s feedback. Students were given feedback on the quality of their
programs based on the model. They were given the chance to correct their
programs after receiving feedback and have it manually reassessed.

###  3.1 Dataset

The dataset was collected over three years from an introduction to computer
science course which teaches Python 3. It consists of four separate
programming assignments which involve a wide range of programming skills. The
simplest is “Travel,” an assignment that involves the distance a vehicle
travelled after going a constant speed for a specified duration. There are 118
student programs for “Travel.” The next assignment, “Budget” is a budgeting
program that lets a user specify a budget and expenses and determines if they
are over or under their budget. 168 student programs were collected for
“Budget”. The third assignment in the dataset consists of creating a program
to play “Rock-Paper-Scissors” against the computer. For this assignment, there
are 111 student programs. The last assignment is programming the classic
casino game “Craps” which involves rolling multiple dice and placing different
types of bets and wagers. This assignment has 120 collected student programs.

All the assignments require the student to write the program from scratch in
Python 3. The programs are to have a command-line, text-based interface and
user validation. Students are required to use if-statements, loops, user
defined functions. The “Craps” program also requires the student to do
exception handling, and file I/O. A requirement of the program was to maintain
a record of their winnings across sessions of playing the game, hence the
results were required to be stored to a file. Also, the standards of design
quality go up as the course progresses and since “Craps” is the last
assignment, it has the highest standards. Each student program has an
associated design score that was normalized to value between zero and one.

###  3.2 Baseline Methods

The model is compared against a variety of baseline regression methods. The
simplest is linear regression, which simply learns a weight per each feature.
Next is a regression decision tree which is trained with the CART algorithm
[1]. It has the advantage over linear regression in that it can learn non-
linear relationships. Non-linearity means a model can learn “sweet-spots”
rather than simply having a “more is better” understanding of some features.
For example, having some modularity in the form of user defined functions is
good, however, too many is cumbersome. The “correct” number of user defined
functions likely should fall into a relatively small range. Model selection on
the maximum depth of the tree with a development set was used to determine
that 10 was the best setting.

However, both of these models have the issue that they are not constrained to
produce a score between zero and one, their prediction can be any real number.
Hence another baseline method was used, created to be an intermediary step
between linear regression and the neural network model. It is a linear model
with a sigmoid transformation which guarantees the output be between zero and
one. This model is effectively the final layer of the neural network model,
i.e. the neural network without the hidden layer. The model is specified by
the equation:

$\displaystyle y=\sigma\bigg{(}\sum\limits_{i=1}^{d}w_{i,j}x_{i}\bigg{)}$ (5)

This model is also trained with ADAM [7]. All the baseline models and the
neural network are trained with MSE as the loss function.

###  3.3 Results

The model was compared versus each baseline in five different settings: the
“Travel”, “Budget”, “Rock-Paper-Scissors”, and “Craps” programs, and a
combined dataset which includes all the programs. The results of the
experiments can be found in Table 1. Each model is evaluated according to two
different metrics: MSE and average accuracy. Average accuracy is defined as
$\frac{1}{n}\sum\limits_{k=1}^{n}\big{(}1-|y_{k}-y^{*}_{k}|\big{)}$.

Overall, the decision tree and the sigmoid-transformed were clearly the two
worst models. This was surprising since decision trees are generally thought
to be strictly more powerful than linear models. However, decision trees look
for highly discriminative features to partition the data into more consistent
groups. The under-performance of the decision tree possibly indicates that
none of the features were especially indicative of a good or bad design on
their own. Instead, the quality of a program is better described by a
collection of subtle features, which gives credence to the belief that design
quality is nuanced.

The reason sigmoid-transformed linear model under-performed linear regression
was likely due to it being trained with ADAM. ADAM does not guarantee
convergence to a global optimum like the analytical solution to linear
regression. Apparently the restriction on predictions to be within the
specified range of zero to one was not important.

Linear regression did surprisingly well, beating both the neural network and
network ensemble in the “Craps” and combined datasets, though barely. In those
cases, the difference between the ensemble and linear regression was less than
a percent. This is likely due to the stability of linear regression’s
predictions. Though linear regression does not have the power and flexibility
of neural networks this can also be a benefit by limiting how wrong their
predictions are. Neural networks and even ensembles can make overconfident
predictions on outliers or other unusual cases. The “Craps” dataset contained
the most complex programs and it is likely a handful of predictions
significantly brought down the average.

The network ensemble outperformed the single neural network in every case,
which is to be expected. The margin of improvement of the ensemble versus the
single neural network in accuracy on the four individual program datasets
ranged from 1% to 4%. The network ensemble did the best overall by being the
best in most cases or coming in a close second in all the other cases. The
importance of using an ensemble is evident on the “Craps” dataset where the
individual neural network under-performed significantly. On the other
datasets, the neural network outperformed linear regression by a small margin,
but on “Craps” the neural network model under-performed the linear regression
model by 5%. Again, this is most likely due the instability and variability of
neural network predictions i.e. small differences in features can lead to a
large difference in the prediction. In the “Craps” dataset, the improvement of
the ensemble over the single neural network model illustrates the relative
stability of the ensemble’s predictions. In every case, the ensemble is
superior to the single neural network and had the best overall performance by
producing the most accurate results on three of the datasets and effectively
tying for the best on the other two.

One noticeable pattern was that all the models performed better on the
“Travel” and “Budget” datasets than on the “RPS”, “Craps”, and combined
datasets. Universally, the most difficult dataset was “Craps” which likely
lowers the accuracy on the combined dataset. Due to the shifting standards and
expectations of student assignments, a model per assignment appears to
worthwhile. This is a bit counter-intuitive since there are many common
standards and expectations across assignments.

Overall, the network ensemble produced reliable, accurate results when trained
per dataset. The accuracy of the ensemble is arguably close to being useful in
practical application. Further, comparing the scores of an instructor versus
another instructor or even against themselves, the rate of agreement must be
less than 100% and with an accuracy of the network ensemble ranging from
83.67% to 94.27%, the model’s accuracy is possibly close to a realistic
ceiling.

###  3.4 Feedback Study

In order the evaluate the effectiveness of the personalized feedback, we
conducted a small study on the effect of the feedback on the design score of
student programs. For the “Rock-Paper-Scissors” program the network ensemble
was used to generate personalized feedback for the student programs instead of
the usual instructor feedback. The network ensemble was the same as used in
the design score experiments, it was trained with prior years worth of student
programs. Having received the personalized feedback, students opted into
correcting their program for extra credit on their assignment. The feedback
was in form of a series of comments, where each comment was “increase” or
“decrease” the name of a feature as described in Section 2.4.

The class is an introduction to computer science course with multiple sections
and two different instructors. Students from both instructors participated in
the study. Out of 73 students enrolled across the sections of the course, 15
students chose to opt-in.

The revised programs were assessed again manually for design quality and the
scores were compared against the originals. The design score of the programs
started at an average of 68.33% and after the feedback and correction the
average rose to 87.92%, a 19.58% absolute improvement. Using a paired t-test,
the improvement was judged to be significant with a p-value of 0.001.

The results of the study suggest the feedback was generally useful in guiding
students to improve the design quality of their programs. The improvement was
noticeable to the instructors anecdotally as well. For example, the usage of
global variables and “magic numbers” decreased significantly. Though the study
does have some caveats including its small sample size and opt-in
participation. It could be that those students willing to opt-in are those
most willing or able improve with a second chance.

##  4 Conclusions & Future Work

Overall, we proposed a neural network model ensemble for predict the design
quality score of a student program and experimentally demonstrated its
effectiveness. Further, our system provided personalized feedback based on the
difference between a program’s feature values and the average features’ value
of “good” programs. A small study provides evidence that the feedback was of
practical use to students. Students were able to improve their programs
significantly based on the feedback they received.

There is also evidence that training models per assignment is most effective.
However, the model needs to be evaluated on more programming assignments.
Further, there is a possibility of utilizing transfer learning [10] to help
the model learn what is in common across the assignments.

The feedback given was shown to be effective, but more nuanced feedback could
be useful. Specifically, feedback targeted to individual lines or segments of
code would possibly help students improve their program’s more effectively.
However, this may require additional supervision i.e. annotation for explicit
training. Active learning [14] or multi-instance learning [3] may be
alternatives to gathering additional annotation.

## References

  * [1] L. Breiman, J. Friedman, C. J. Stone, and R. A. Olshen.  Classification and regression trees.  CRC press, 1984.
  * [2] Code.org.  Code.org: Learn computer science.  https://code.org/research.
  * [3] T. G. Dietterich, R. H. Lathrop, and T. Lozano-Pérez.  Solving the multiple instance problem with axis-parallel rectangles.  Artificial intelligence, 89(1-2):31–71, 1997.
  * [4] A. Efremov, A. Ghosh, and A. Singla.  Zero-shot learning of hint policy via reinforcement learning and program synthesis.  International Educational Data Mining Society, 2020.
  * [5] A. Gupta and N. Sundaresan.  Intelligent code reviews using deep learning.  In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD’18) Deep Learning Day, 2018.
  * [6] G. E. Hinton, N. Srivastava, A. Krizhevsky, I. Sutskever, and R. R. Salakhutdinov.  Improving neural networks by preventing co-adaptation of feature detectors.  arXiv preprint arXiv:1207.0580, 2012.
  * [7] D. Kingma and J. Ba.  Adam: A method for stochastic optimization.  In the 3rd International Conference on Learning Representations (ICLR), 2014.
  * [8] V. Nair and G. E. Hinton.  Rectified linear units improve restricted boltzmann machines.  In Proceedings of the 27th international conference on machine learning (ICML-10), pages 807–814, 2010.
  * [9] B. Paassen, B. Hammer, T. W. Price, T. Barnes, S. Gross, and N. Pinkwart.  The continuous hint factory-providing hints in vast and sparsely populated edit distance spaces.  Journal of Educational Data Mining, 2018.
  * [10] S. J. Pan and Q. Yang.  A survey on transfer learning.  IEEE Transactions on knowledge and data engineering, 22(10):1345–1359, 2009.
  * [11] C. Piech, J. Huang, A. Nguyen, M. Phulsuksombati, M. Sahami, and L. Guibas.  Learning program embeddings to propagate feedback on student code.  In International conference on machine Learning, pages 1093–1102. PMLR, 2015.
  * [12] Pylint.org.  Pylint - code analysis for python.  https://pylint.org.
  * [13] F. Rosenblatt.  Principles of neurodynamics. perceptrons and the theory of brain mechanisms.  Technical report, DTIC Document, 1961.
  * [14] B. Settles.  Active learning literature survey.  Technical report, University of Wisconsin-Madison Department of Computer Sciences, 2009.
  * [15] G. Van Rossum, B. Warsaw, and N. Coghlan.  Pep 8.  https://www.python.org/dev/peps/pep-0008/.
  * [16] L. Wang, A. Sy, L. Liu, and C. Piech.  Learning to represent student knowledge on programming exercises using deep learning.  International Educational Data Mining Society, 2017.

##  Appendix A Feature Functions

  * •

The number of functions

  * •

The number of assignments

  * •

AST nodes per function

  * •

Lines of code per function

  * •

Total lines of code

  * •

Number of literals

  * •

The proportion of white-space characters to the total number of characters

  * •

Number of empty lines

  * •

Deepest level of indentation

  * •

Number of “if” statements

  * •

Number of comments

  * •

Number of AST nodes per lines of code

  * •

Number of try-except statements

  * •

AST nodes per try-except statement

  * •

AST nodes per “if” statement

  * •

Number of lists

  * •

Number of tuples

  * •

Average line number of literals

  * •

Average line number of function definition

  * •

Average line number of “if” statement

  * •

Ratio of AST nodes inside functions versus total number of AST nodes

  * •

Number of function calls

  * •

Number of “pass” statements

  * •

Number of “break” statements

  * •

Number of “continue” statements

  * •

Number of global variables

  * •

Number of zero and one integer literals

  * •

Average line number of “import” statement

  * •

Number of numeric literals

  * •

Number of comparisons

  * •

Number of “return” statements

  * •

Maximum number of “return” statements per function

  * •

Maximum number of literals per “if” statement


# kafe2 - a Modern Tool for Model Fitting in Physics Lab Courses

Johannes Gäßler(1), Günter Quast(1), Daniel Savoiu(2), Cedric Verstege(1)
(1) Karlsruhe Institute of Technology (KIT), (2) Hamburg University

###### Abstract

Fitting models to measured data is one of the standard tasks in the natural
sciences, typically addressed early on in physics education in the context of
laboratory courses, in which statistical methods play a central role in
analysing and interpreting experimental results. The increased emphasis placed
on such methods in modern school curricula, together with the availability of
powerful free and open-source software tools geared towards scientific data
analysis, form an excellent premise for the development of new teaching
concepts for these methods at the university level. In this article, we
present kafe2, a new tool developed at the Faculty of Physics at the Karlsruhe
Institute of Technology, which has been used in physics laboratory courses for
several years. Written in the Python programming language and making extensive
use of established numerical and optimization libraries, kafe2 provides simple
but powerful interfaces for numerically fitting model functions to data. The
tools provided allow for fine-grained control over many aspects of the fitting
procedure, including the specification of the input data and of arbitrarily
complex model functions, the construction of complex uncertainty models, and
the visualization of the resulting confidence intervals of the model
parameters.

##  1 Model fitting in lab courses: pitfalls and limitations

A primary goal of data analysis in physics is to check whether an assumed
model adequately describes the experimental data within the scope of their
uncertainties. A regression procedure is used to infer the causal
relationships between independent and dependent variables, the abscissa values
$x_{i}$ and the ordinate values $y_{i}$. The dependency is typically
represented by a functional relation $y_{i}=f(x_{i};\\{p_{k}\\})$ with a set
of model parameters $\\{p_{k}\\}$. Both the $y_{i}$ and the $x_{i}$ are the
results of measurements and hence are affected by unavoidable measurement
uncertainties that must be taken into account when determining confidence
intervals for the model parameters. It has become a common standard in the
natural sciences to determine confidence intervals for the parameter values
based on a maximum-likelihood estimation (MLE). In practice, one typically
uses the negative logarithm of the likelihood function (NLL), which is
minimised by means of numerical methods.

In the basic laboratory courses, parameter estimation is traditionally based
on the method of least squares (LSQ), which in many cases is equivalent to the
NLL method. For example, fitting a straight line to data affected by static
Gaussian uncertainties on the ordinate values $y_{i}$ is an analytically
solvable problem and the necessary calculations can be carried out by means of
a pocket calculator or a simple spreadsheet application. For simplicity, non-
linear problems are often converted to linear ones by transforming the $y$
values. However, the Gaussian nature of the uncertainties is lost in this
process, which is why the transformed problem only gives accurate results if
the uncertainties are sufficiently small. Simple error propagation is often
used to transform the fitted parameters like the slope or $y$ intercept to the
model parameters of interest. As a result, the parameter uncertainties
obtained in this way are no longer Gaussian, and it is non-trivial to
determine sensible confidence intervals. When the independent variables are
also affected by measurement uncertainties, or when relative uncertainties
with respect to the true values are present, this traditional approach reaches
its limits.

Such simple analytical methods are very limited in their general applicability
to real-world problems. For instance, the presence of uncertainties on the
abscissa values in addition to those on the ordinate values, as they would be
encountered when measuring a simple current-voltage characteristic of an
electronic component, already results in a problem that can generally no
longer be solved analytically. Unfortunately, most of the widely used
numerical tools for parameter estimation provide only limited support for such
non-linear problems, even though they occur frequently in practice.

Prior to drawing any conclusions on model parameters, the most important task
in physics is to check the validity of a model hypothesis assumed to describe
a set of measurement data. The full task therefore consists of the following
steps:

  1. 1.

carefully quantifying the uncertainties of all input measurements;

  2. 2.

defining a model hypothesis for the measured data;

  3. 3.

testing whether the measurements are compatible with the model hypothesis;

  4. 4.

if so, determining the model parameters, e.g. the slope and intercept of a
straight line.

It is precisely the hypothesis testing in step 3 that drives progress in
scientific understanding by providing the relation between measurements and
theoretical models. A simple example for such a hypothesis test is the
$\chi^{2}$ test, which results directly from the frequently used method of
least squares. A test for the validity of the model assumption should
definitely be performed before drawing any conclusions from the values and
uncertainties of the fitted model parameters. Unfortunately, many of the
common tools come with a default setting that assumes the correctness of the
given parameterization and scales the parameter uncertainties so that the
model perfectly describes the data within the parameter uncertainties so
determined. This behaviour can usually be switched off by choosing appropriate
options, but in practice this is often neglected.

The parameters of a complex model can be highly correlated. In this case it is
no longer sufficient to consider the distribution of just a single parameter.
Instead the shared distribution of multiple parameters needs to be considered,
for example when conducting hypothesis testing. Correlations between model
parameters can also be problematic for numerical optimisation. For this reason
strategies for choosing appropriate parameterizations to reduce these
correlations also need to be addressed. It is then indispensable that the
correlations are also shown when presenting the fit results. Ideally the
correlation between two model parameters can be expressed with a simple
correlation coefficient. However, if non-linear models are fitted, this is
often not sufficient for the evaluation of the result; in this case the
corresponding confidence regions should be determined and presented as contour
graphs for pairs of parameters.

The distinction between ”statistical” and ”systematic” uncertainties leads to
many misunderstandings among students. For example, the systematic error of a
measuring instrument becomes a statistical one if several measuring
instruments of the same type are used. A much more suitable approach is to
differentiate whether uncertainties affect all measured values or groups of
measured values equally or whether they are independent. This approach,
however, requires dealing with the covariance matrix of measurement
uncertainties and constructing it from the individual uncertainties on a
problem-by-problem basis. Unfortunately, there are very few simple tools that
allow a full covariance matrix to be considered in the fitting process. None
of the tools commonly used in physics lab courses allow students to take into
account correlated uncertainties of both the abscissa and ordinate values.

While there are tools with which the requirements listed here could be
partially implemented, these do not provide the simplicity needed for
undergraduate physics education. To address the issues described above it
became necessary to write a tool of our own, kafe2, an open-source Python
package designed to provide a flexible Python interface for the estimation of
model parameters from measured data.

##  2 The Open-Source Python Package kafe2

After designing a prototype package, kafe, development continued with kafe2
[1]. Said tools make use of contemporary methods for visualization and
evaluation of measurement data and are based on freely available, open-source
software that is also used in scientific contexts and thus relevant for later
professional practice.

###  2.1 Implementation

Due to its widespread use in scientific communities as well as in the
burgeoning professional field of data science, _Python_ [2] was chosen as the
programming language for kafe2. As a high-level language with particular
emphasis on clear and intuitive syntax, and convenient features such as
dynamic typing and garbage collection, Python provides a beginner-friendly
programming environment, requiring comparatively little knowledge to use
effectively. This is particularly advantageous in a field such as physics,
which requires a very broad foundation of knowledge prior to specialization
and thus allows only a limited percentage of the curriculum to be assigned to
programming.

However, an inevitable downside of Python is that it is comparatively slow
compared to low-level languages such as C/C++. Fortunately, the Python
ecosystem features many libraries for numerical computation, most notably the
_NumPy_ [3] library, which makes use of fast underlying C implementations for
running such computations efficiently. Optimal use of these libraries does
require some detailed knowledge, but in the case of kafe2 this generally only
applies to development rather than use.

The numerical optimization algorithms necessary for finding the global minimum
of the negative log-likelihood function are provided either by the library
_SciPy_ [5] (”scientific Python”) or by the Python package _iminuit_ [6],
which is based on the _MINUIT_ package developed at CERN, the European
Organization for Nuclear Research.

The graphics library _matplotlib_ [4] is used as the graphical backend for
publication quality visualizations of input data and results. kafe2 itself is
implemented in pure Python and can therefore run on all common platforms with
a Python environment, including ARM-based platforms such as the Raspberry Pi.

###  2.2 Special Features of kafe2

The kafe2 package provides numerous unique features that are listed in the
following.

To represent complex uncertainty models with absolute and/or relative,
independent and/or correlated uncertainties on abscissa and/or ordinate
values, kafe2 uses two covariance matrices collecting the uncertainty
components on the abscissa and ordinate directions. As described below, these
are combined into a single overall covariance matrix used in the numerical
regression. Users can specify multiple uncertainties of the same kind which
are then automatically combined into the corresponding covariance matrices. If
relative uncertainties are specified, they most often refer to the true values
of the input data, which are only (approximately) known after the regression.
kafe2 is able to handle such cases by dynamically calculating the
uncertainties as relative to the model; optionally and for compatibility with
other tools it is possible to use relative uncertainties with respect to the
measured data. Correlations between uncertainty components or between
individual uncertainties can be specified as a correlation coefficient or as a
complete covariance matrix.

In addition to pre-defined, rather simple model functions, arbitrarily complex
user-defined functions can be easily implemented using native Python syntax.
It is also possible to simultaneously fit several model functions with shared
parameters to data; this feature is very useful if auxiliary or external
measurements of model parameters need to be considered. Fixing the values of
one or more model parameters or constraining them within given external
uncertainties is also supported.

The numerical fitting procedure in kafe2 relies on the NLL method, but falls
back to the simple LSQ method if equivalent to NLL in order to avoid
computationally expensive operations. Distributions can be fitted to binned
data (histograms) using a Poisson likelihood; optionally and for compatibility
with other packages the LSQ method can alternatively be chosen. kafe2 also
supports fits with user-defined maximum-likelihood functions.

Parameter uncertainties are determined in the traditional way relying on the
the Cramér-Rao-Fréchet bound, i.e. from the second derivatives of the LSQ or
NLL function at the minimum. To deal with non-linear cases and correlated
parameters, a likelihood-based method known as ”profile likelihood” is also
optionally offered for determining parameter confidence intervals. In multi-
parameter problems, this is also useful to account for the influence of
nuisance parameters when determining individual parameter uncertainties.

Fit results are provided in graphical form, but output in text form is also
possible. Access functions exist to extract results for further processing
with Python code. If meaningful, kafe2 provides metrics to judge the goodness
of a fit result. The output also includes the parameter correlations, as is
standard for most other fitting tools. The graphical output shows the data and
the best-fit model function with an optional representation of the model
uncertainty as a shaded band. An optional graphical representation of the
profile likelihood for each parameter and of confidence regions of pairs of
parameters as contours are also provided.

More details on the mathematical procedures and a concrete application example
are given below.

###  2.3 Mathematical procedures

The estimation of model parameters and their uncertainties in kafe2 is based
on the method of maximum-likelihood estimation. Given a set of $N$
measurements $y_{i}$ with corresponding abscissa values $x_{i}$, as well as a
model function $f({x_{i}},\\{p_{k}\\})$ depending on a set of parameters
$\\{p_{k}\\}$, the most commonly assumed data distribution is the
$N$-dimensional Gaussian distribution $\cal N$ with the uncertainties on the
ordinate values being described by the covariance matrix $\bf V$:

$\mathcal{N}\left(\\{y_{i}\\},f(\\{x_{i}\\},\\{p_{k}\\}),{\bf
V}\right)\,=\,\displaystyle\frac{1}{\sqrt{(2\pi)^{N}\det{\bf
V}}}\exp\left(-\frac{1}{2}\sum_{i,j=1}^{N}\left(y_{i}-f(x_{i},\\{p_{k}\\})\right)\left({\bf
V}^{-1}\right)_{ij}\left(y_{j}-f(x_{j},\\{p_{k}\\})\right)\right)\,.$

The cost function ${\cal{C}}$ to be minimized is defined as twice the negative
logarithm of the likelihood ${\cal{L}}$ of $\cal N$,
${\cal{C}}=-2\ln{\cal{L}}$. Terms that are independent of the parameters are
typically neglected in the definition of the likelihood since they are
irrelevant for the location and shape of the global cost function minimum.

The total uncertainties of the measurement data are described by the
covariance matrix ${\bf V}$. The uncertainties of the measurement points in
the abscissa direction are described by a covariance matrix with elements
${\bf V}^{x}_{ij}$. Because the Gaussian likelihood is only defined for
ordinate uncertainties kafe2 handles abscissa uncertainties by transforming
them into ordinate uncertainties. This is achieved by multiplying them with
the first derivatives of the model function. The resulting uncertainties are
then simply added to the covariance matrix elements ${\bf V}^{y}_{ij}$ of the
ordinate uncertainties. The elements of the total covariance matrix used in
the fit are thus

${\bf V}_{ij}={\bf V}^{y}_{ij}+{\bf V}^{x}_{ij}\cdot f^{\prime}(x_{i})\cdot
f^{\prime}(x_{j})\,.$

This approach assumes that, near the global cost function minimum, the fitted
model function can be sufficiently well approximated by a straight line in the
vicinity of the support points $x_{i}$, i.e. that second-order derivatives are
negligible on the scale of the abscissa uncertainties $\Delta x_{i}$.

For uncertainties that are independent of the parameter values, ${\cal{C}}$ is
equivalent to the method of least squares, which is classically used as the
cost function in most fitting algorithms,

${\cal{C}}(\\{y_{i}\\},{\bf
V}^{y},f(\\{x_{i}\\},\\{p_{k}\\}))\,=\,\displaystyle\sum_{i,j}\left(y_{i}-f(x_{i},\\{p_{k}\\})\right)\left(\big{(}{{\bf
V}^{y}}\big{)}^{-1}\right)_{ij}\left(y_{j}-f(x_{j},\\{p_{k}\\})\right)\,.$

The expected value of ${\cal{C}}$ for the best parameter values estimated from
a hypothetical set of data follows a $\chi^{2}$-distribution with a number of
degrees of freedom (NDF) corresponding to the difference between the number of
data points and the number of free model parameters. The ratio
$\mathcal{C}/\mathrm{NDF}$, often called the ”reduced $\chi^{2}$”, has an
expectation value of unity for a model that accurately describes the data
within their uncertainties. The ”$\chi^{2}$ probability” $p=1-q$ is determined
from the quantile $q$ of the $\chi^{2}$ distribution for the observed minimum
and can be used to quantify the goodness-of-fit. If, however, the
uncertainties do depend on the parameters, as is for instance the case for
relative ordinate uncertainties or uncertainties on the abscissa values, then
the method of least squares is no longer equivalent to the method of maximum
likelihood. This is because the normalization term of the multivariate
Gaussian, which is a function of $\det\left(\bf V(\\{p_{k}\\})\right)$, must
also be taken into account. While a parameter-dependent covariance matrix
needs to be recomputed on every parameter change, thus adding to the
computational cost of the minimization, it remains feasible on modern
hardware, particularly if the inversion of the covariance matrix $\bf V$ is
replaced by a Cholesky decomposition ${\bf L}{\bf L}^{T}={\bf V}$. In
addition, the Cholesky decomposition also enables the quick computation of the
covariance matrix determinant by exploiting that $\bf L$ is a triangular
matrix: $\det\left({\bf V}\right)=\det\left({\bf L}\right)\cdot\det\left({\bf
L}^{T}\right)=2\det\left({\bf L}\right)$.

Following common standards based on the Cramér-Rao-Fréchet bound, an estimate
of the covariance matrix of the parameter uncertainties is determined from the
inverse Hessian matrix, containing the second-order derivatives of
${\cal{C}}(\\{y_{i}\\},V,f(\\{x_{i}\\};\\{p_{k}\\}))$ with respect to the
parameters $p_{k}$. In general, however, this approach is not sufficient for
reliably determining the confidence intervals of the fitted parameters. For
this reason, kafe2 also implements the _profile likelihood_ method. The so-
called profile as a function of a subset of the parameters is determined by
keeping the values of these parameters fixed and minimizing ${\cal{C}}$ with
respect to all other parameters. The boundaries of a confidence interval are
determined by finding the locus of points for which the profile likelihood
increases by a value $\Delta\mathcal{C}$ with respect to the global minimum.
In the case of two parameters, this amounts to finding the intersection of the
two-dimensional profile likelihood function with a horizontal plane above the
global cost function minimum.

The value of $\Delta{\cal{C}}$ determines the confidence level of the
confidence interval/region. For example, for a one-dimensional confidence
interval $\Delta{\cal{C}}=1$ corresponds to the $1\,\sigma$ interval of a
Gaussian distribution with a confidence level of 68.3%, where $\sigma$ is the
standard deviation of the one-dimensional Gaussian distribution. In the two-
dimensional case (as in the profile likelihood of two parameters) the
intersection of a two-dimensional profile likelihood with a horizontal plane
is a one-dimensional curve. This curve corresponds to a confidence contour,
which is an ellipse in case of Gaussian uncertainties111 It should be noted
that the uncertainties on the fitted model parameters can often be made more
Gaussian-like by changing the parameterization of the model function.. If the
contour deviates strongly from the elliptical form, is no longer possible to
correctly specify a confidence region of the two parameters by providing a
simple correlation coefficient. In this case, the contour plots of the
affected pairs of parameters as well as the one-dimensional confidence
intervals should be documented.

The profile likelihood method is important because non-linear problems are
very common in practice. Even a simple linear regression with additional
uncertainties on the abscissa values leads to a non-linear problem for which
no general analytical solution exists. Relative uncertainties on the ordinate
values introduce similar issues. They should be specified relative to the
model (as a stand-in for the unknown ”true values”) to avoid the bias which
otherwise would give too much weight to data points with downward fluctuations
of the measured values and hence the assigned uncertainties. When the model
function itself is nonlinear in the parameters or when abscissa uncertainties
or relative ordinate uncertainties are used, the shape of the covariance
contours should always be inspected. Using appropriate tools, the associated
calculations usually only take on the order of seconds on consumer hardware.
An example with kafe2 is shown below.

In addition to the Gaussian-based cost functions discussed until now, kafe2
also offers other cost functions, for example a cost function based on the
Poisson distribution for fitting models to frequency distributions
(histograms). The specification of user-defined cost functions is also
possible; in order to benefit from the statistical interpretation of the fit
results provided by kafe2, such cost functions must be based on valid
likelihood functions.

Often it is desirable to show the uncertainty of the fitted model function. An
approximate confidence region $\Delta f(x;p)$ can be derived from the
covariance matrix, ${\bf V}^{p}$, of the model parameters and the gradient,
$\nabla_{p}f(\\{x_{i}\\};\\{p_{k}\\})$, of the model function with respect to
the parameters by (linear) propagation of the uncertainties:

$\Delta
f(x_{i},\\{p_{k}\\})=\Big{(}\nabla_{p}f(x_{i};\\{p_{k}\\})\Big{)}^{T}{\bf
V}^{p}\ \nabla_{p}f(x_{i};\\{p_{k}\\})\,.$

##  3 User interface and applications

To mitigate the inherent trade-off between usability and complex functionality
kafe2 offers the user several interfaces for fitting:

  * •

For those without programming knowledge, or for the convenience of advanced
users, a command line interface “ _kafe2go_ ” is provided. Users only need to
supply a configuration file in the format of the widely used data
serialization language YAML [8], which specifies the data, the model function
and fitting options.

  * •

More flexibility is offered to slightly advanced users with basic programming
knowledge; they can call the Python interface of built-in kafe2 pipelines as
part of a larger Python script.

  * •

Advanced users have the possibility to construct custom pipelines by directly
instantiating the kafe2 objects that represent input data, model functions,
fits or graphical output.

As an example, the stand-alone application kafe2go can be used via the
following command line call:

`kafe2go <name>.yaml`.

Figure 1 illustrates the result of fitting a straight line to data points with
uncertainties in both the ordinate and abscissa directions. The input file
contains only the following lines:

     # input data straight line fitting
     x_data: [1.0, 2.0, 3.0, 4.0]
     x_errors: 5%
     y_data: [2.3, 4.2, 7.5, 9.4]
     y_errors: 0.4

Note: if no model function is specified, kafe2go defaults to a first-degree
polynomial as the model.

Figure 1: Graphical output of kafe2. The shaded region around the fitted
straight line shows the 68% confidence interval of the model function.

When using kafe2 as a Python library, a Jupyter [9] environment is very
convenient. One advantage is that this allows users to combine text, program
code and outputs in a so-called ”Jupyter notebook”. A second advantage is that
it also enables teaching staff to provide the entire Python environment for
running kafe2 as a service to students; the students themselves only need to
connect to a Jupyter server via web browser without the need to spend time
setting up a Python environment themselves.

To use the built-in kafe2 pipelines for model fitting, users simply need to
call the corresponding functions. The creation of custom pipelines, however,
requires a more advanced understanding of object-oriented programming. First,
a data container object is instantiated from numerical input data and their
uncertainties. Then a fit object is created from a data container object and a
model function. After possibly constraining or limiting parameters via
dedicated method calls, the method `do_fit()` is called to perform the
numerical optimization of the cost function. The fit results are then accessed
via the corresponding attributes of the fit object or printed to console by
calling the `report()` method. In order to create a graphical representation
of the fit results a plot object can be instantiated from the fit object. A
complete code example, which corresponds exactly to the one just discussed and
leads to the same graphical output as in Figure 1, looks as follows:

        # Python code for fitting a straight line
        from kafe2 import XYContainer, Fit, Plot
        xy_data = XYContainer(
            x_data=[1.0, 2.0, 3.0, 4.0],
            y_data=[2.3, 4.2, 7.5, 9.4]
        )
        xy_data.add_error(
            axis=’x’,
            err_val=0.05,
            relative=True
        )
        xy_data.add_error(
            axis=’y’,
            err_val=0.4
        )
        line_fit = Fit(xy_data)
        line_fit.do_fit()
        line_fit.report()
        plot = Plot(fit_objects=line_fit)
        plot.plot()
        plot.show()

Since fits of models to two-dimensional data points is a very common use case
a built-in function implementing the pipeline is also available:

        import kafe2
        kafe2.xy_fit(
            x_data=[1.0, 2.0, 3.0, 4.0],
            y_data=[2.3, 4.2, 7.5, 9.4],
            x_error_rel=0.05,
            y_error=0.4
        )
        kafe2.plot()

Because relative abscissa uncertainties are specified (via `x_error_rel`), the
built-in pipeline automatically switches to calculating parameter
uncertainties with the profile likelihood. A plot of the confidence
intervals/regions is also created automatically.

Users are provided with several resources to familiarize themselves with the
use of kafe2. In addition to the traditional documentation explaining the use
by showing examples there is also a tutorial in the form of a Jupyter notebook
that takes a more interactive approach. Among others, the following topics are
covered:

  1. (a)

Simple linear regression with a straight line as model function;

  2. (b)

definition of arbitrary model functions;

  3. (c)

specification of uncertainties and their correlations for both ordinate and
abscissa values;

  4. (d)

treatment of relative uncertainties;

  5. (e)

incremental construction of covariance matrices from individual uncertainties;

  6. (f)

calculation of confidence intervals and contours using the profile likelihood
method;

  7. (g)

constraining model parameters in accordance with external information;

  8. (h)

fitting distributions to binned or unbinned data using the maximum likelihood
method;

  9. (i)

fitting models to data with non-Gaussian uncertainties, e.g. frequency
distributions;

  10. (j)

simultaneous fitting of multiple models with shared parameters to several data
sets.

Since kafe2 is open-source software, the code can be modified or extended as
needed.

###  3.1 Application example

Figure 2 shows the results obtained when fitting both a linear and an
exponential model to the same data. Both the ordinate and the abscissa values
are subject to measurement uncertainty. Additionally the uncertainties on the
ordinate values are specified as relative. Because of the reference to the
model value, the uncertainties shown in the graph differ depending on the
model function.

When plotting or printing fit results, kafe2 provides the user with two
metrics that can be used for hypothesis testing, specifically for Pearson’s
$\chi^{2}$ test. The first is the value of the cost function at the minimum
divided by the number of degrees of freedom $\mathrm{NDF}$, which has an
expectation value of unity, assuming that the model adequately describes the
data within its statistical fluctuations. The second metric quantifying the
goodness-of-fit is the $\chi^{2}$ probability, i.e. the probability to find a
minimum $\chi^{2}$ value larger than the one actually observed, with larger
values indicating a better agreement between the data and the model. Based on
this value, a simple hypothesis test can be performed, rejecting the
hypothesis that the data is adequately described by the model if the
$\chi^{2}$ probability is less than a pre-defined significance threshold
$\alpha$.

In the case shown here, the probability for obtaining an even higher than
observed $\chi^{2}$ cost function value is 8.8% for the linear model and 44.5%
for the exponential model, indicating that both models are acceptable with a
significance of $\alpha=5\%$. As the visual impression already shows, the
exponential model fits the data more closely.

Figure 2: Fitting two models to the same data with kafe2.

Figure 3: Two-dimensional confidence contours and profile likelihood curves
for the exponential model.

The parameterization of the exponential model is - intentionally - chosen
suboptimally to amplify the nonlinearity. The parameter $x_{0}$ is in the
denominator of the expression in the exponent; if $x_{0}$ were replaced by a
parameter $\lambda:=1/x_{0}$, then the parameter uncertainties derived from
the Cramér-Rao-Fréchet bound would be more accurate. The deviation from the
parabolic shape can be inspected via the profile likelihood curve, as shown in
Figure 3. The profile likelihood for the parameter $x_{0}$ \- drawn as a solid
blue line - clearly deviates from the parabola that is extrapolated from the
second derivative at the cost function minimum. The 68.3% confidence interval
is determined by finding the intersection of the likelihood curves with a
constant function equal to 1; the confidence interval derived from the profile
likelihood method differs from the conventional confidence interval obtained
by extrapolating the second derivatives. Instead of a symmetric uncertainty of
$\pm$0.70, this results in an asymmetric confidence interval from -0.66 to
+0.76 around the central value.

The exponential model is a linear function of the parameter $A_{0}$. Despite
this, the profile likelihood of $A_{0}$ also differs from a parabola: the
68.3% confidence interval ranges from -0.27 to +0.32 relative to the central
value. This effect is mostly caused by the relative uncertainties on the
ordinate values, which are specified as relative to the model values. As such,
larger model values increase the uncertainties and therefore result in a lower
cost function value. Conversely, smaller model values result in a higher cost
function value. The asymmetry of the uncertainty of parameter $a$ shown as the
result of the linear fit in Figure 2 is caused by the same effect.

The third graph in the bottom left corner of Figure 3 shows the boundaries of
the confidence regions for the shared distribution of $x_{0}$ and $A_{0}$, as
determined by the profile likelihood method. The plot shows the confidence
contours for $\Delta\chi^{2}=1$ and $\Delta\chi^{2}=4$, which correspond to
the $1\,\sigma$ and $2\,\sigma$ uncertainty ranges, respectively. The slope of
the contour reflects the strong correlation of 90% between the parameters
$A_{0}$ and $x_{0}$. The cross indicates the uncertainties determined from the
second derivatives at the minimum. The deviation from the elliptical shape is
visible for both contours, but it is particularly pronounced for the
$2\,\sigma$ contour. This is expected since the extrapolation of the second
derivatives is more accurate close to the cost function minimum. Students
should be encouraged to include the asymmetric uncertainties as well as the
confidence contours when documenting their results.

##  4 Experiences and conclusion

The methods upon which kafe2 is built are a core part of the physics
curriculum at the Karlsruhe Institute of Technology. Bachelor students are
required to attend a lecture on computer-based data analysis which includes
topics such as the basics of statistics and the visualization of data. The
lecture is aimed at students in the second semester and serves as a
preparation for the physics laboratory courses starting in the third semester.
As such, the focus of the lecture is placed on equipping students with the
tools and knowledge necessary to perform simple but modern data analyses on
their own. More complex problems such as fits with user-defined maximum-
likelihood functions or fits to histogram data arise in the context of
advanced practical courses or work on projects done as part of a Bachelor’s
thesis.

The kafe and kafe2 packages have been in continuous use at KIT in the context
of physics laboratory courses since 2015. For ease of use, sample code is
provided as part of a collection of useful Python functions for storage,
processing and visualization (software package _PhyPraKit_ [10]), which
encapsulates the complex functionality of kafe2 and provides a simple
interface sufficient for most practical needs. Optional tasks offer incentives
to try out more complex functionality, replacing the traditional simple error
calculation with complete uncertainty models with covariance matrix in the
model fitting.

More in-depth lectures on data analysis that build on the foundations laid in
the practical courses are offered to postgraduate students. Positive feedback
from supervisors of Bachelor’s theses indicates that the students experience a
significant gain in competence when it comes to the evaluation of scientific
data as well as its visualization, discussion, and presentation in scientific
contexts. The student surveys for the courses show that the prior knowledge in
the area of digital data processing varies greatly from student to student.
The acceptance of the tools provided, the assessment of their usefulness, and
the positive awareness of the students’ own competence increase over the
course of their studies.

The lecture on computer-based data analysis is deliberately placed early in
the physics curriculum. Students can only experience a gain in knowledge from
laboratory courses if they understand the mechanisms that the experiments are
based on - this holds true both for the underlying physics and the statistical
methods employed. However, the early introduction of software tools also
necessitates that said tools are (comparatively) easy to use. Most physics
students were taught at most basic programming in school and cannot be
expected to learn how to apply the complex tools used in scientific practice.
At the same time the available simple tools (some of which require no
programming knowledge at all) do not use state-of-the-art statistical methods
and thus cannot enlighten students when it is necessary to apply them. To
brigde this gap, kafe2 aims to provide for an easy-to-use tool built on state-
of-the-art statistical methods. While the early introduction of modern
analysis methods still poses some difficulty, most students respond positively
to the challenge and are ultimately able to master it well.

With the tools commonly used in basic physics education for model fitting and
visualization, an analysis as shown in the examples above is practically
impossible to perform. With kafe2, however, students have a sufficiently
simple tool at their disposal to carry out data analysis tasks in accordance
with contemporary scientific standards.

## References

  * [1] Description and source code kafe2,
URLs: https://philfitters.github.io/kafe2/ and
https://github.com/PhiLFitters/kafe2/.

  * [2] Interpreted programming language Python, URL: https://www.python.org/.
  * [3] Scientific computing library NumPy, URL: https://numpy.org/,
Harris, C.R., Millman, K.J., van der Walt, S.J. et al., Array programming with
NumPy, Nature 585, 357–362 (2020), DOI: 10.1038/s41586-020-2649-2.

  * [4] Graphics library matplotlib, URL: https://matplotlib.org/,
J. D. Hunter, Matplotlib: A 2D Graphics Environment, Computing in Science &
Engineering, vol. 9, no. 3, pp. 90-95, (2007), DOI: 10.5281/zenodo.7162185.

  * [5] Software library for mathematics, science and engineering SciPy, URL: https://www.scipy.org/,
Pauli Virtanen et al., SciPy 1.0: Fundamental Algorithms for Scientific
Computing in Python, Nature Methods, 17(3), 261-272 (2020), DOI:
10.1038/s41592-019-0686-2.

  * [6] iminuit, URL: https://github.com/scikit-hep/iminuit/, Hans Dembinski, Piti Ongmongkolkul et al., DOI: 10.5281/zenodo.3949207.
  * [7] Data analysis framework ROOT, URL: https://root.cern/.
  * [8] Data description language YAML, URL: https://yaml.org/.
  * [9] Interactive computing service Jupyter, URL: https://jupyter.org/.
  * [10] Software package PhyPraKit for digital data processing in physics laboratory courses,
URL: https://github.com/GuenterQuast/PhyPraKit/.