Design of objective functions for model calibration

Introduction

Model calibration or parameter estimation is the task of tailoring the free parameters of the model to best explain the known observations. This is vital because making new predictions from the model follows this step. For the predictions to be accurate the free parameters have to be first identified with the known data. An effective passage through this step requires fighting on multiple fronts. First one is the design of the model itself. This seems trivial to say but its often not easy to test all the behaviors of the model and so it is an iterative process in general. If the starting point is to use a first order polynomial relation and it doesn't accurately fit the observations then it might be better to try higher order non-linear relations.

The second front for effectively attending the calibration task is to choose a search strategy. In case of continuous and differentiable functions calculus gives us the closed form solutions with the use of derivatives, in reality we often encounter functions too complicated to calculate derivatives. In such cases we rely on search/heuristic methods ranging from a simple grid search to advanced meta-heuristic methods (Genetic Algorithm for example).

Because an effective search relies on evaluating the fitness values at different positions in the search space, it is important to choose an appropriate objective function which is the third battle front. An objective function maps a position in space to a value (Scalar or a vector). If the goal of the search is to reach a known destination point the objective function can be thought of as a distance measure. In some cases we don't know the destination exactly but understand its characteristics, like the final point maximizes or minimizes some measure. In either cases designing the objective function that resonates with the particular instance is crucial to reach the ideal solution. For example, Consider a cube whose edges are the path connecting the vertices. If a particle is to move from one vertex to another, then picking an Euclidean distance measure is ineffective over a Manhattan distance measure. Since the particle cannot move along the diagonal, the Euclidean distance doesn't give a true measure of distance between points on opposite sides of a diagonal but a Manhattan distance does.

Path taken by a particle moving on cube along its edges from A to G is better expressed by Manhattan distance(dark edges) than a Euclidean distance(dotted red)

In this post I will share two instances where designing a functional form by incorporating the nuances of the case helped me implement an effective optimization.

Current value vs Cumulative value:

In a path finding problem the goal is to find a path between two points in a given map. For example we want to find a route between our home to a restaurant. To be efficient one would like to follow the shortest path or a thrifty one. The problem gets more realistic if the map contains roadblocks, circular paths, sections with difficult terrain to drive etc. To gauge the efficiency of a path chosen it is necessary to know the distance between points. On a map like a well planned city, the Manhattan distance is a good measure for the same reason as with the cube. But in such a setting one would easily run into the problem of circular path i.e returning to visited points and hence get struck in a loop. So an objective function better than the just distance to goal is needed. In this instance that would be to use a cumulative distance. That means instead of just learning the distance to destination from the current location, it is better to keep track of the distance traveled.

In a restricted map a simple measure of distance from current location to destination can lead to an impasse. A cumulative measure will play the long game and consider more options

An accumulated measure doesn't avoid taking a path that seems costly immediately but rather plays the long game. Such an approach of fitness measure is needed because there is no shortest segment path linking the origin and destination and ignoring this will lead to failure of reaching the destination.

Choice of function based on the magnitude of values

In an experiment undertaken to understand the relationship between two variables, multiple observations are noted and model fitting is used to find the best parameters of the model. For this fitting job a distance measure like a square of difference between the corresponding points is used to measure the fitness of a particular set of parameters. A slightly different but analogous situation is a multi-objective case where there are several measurable/observable variables and the model has to be optimized simultaneously on all measures. This array of measures can be converted to the usual case of scalar output by finding appropriate weights. These weights deserve special attention when combining measures with different ranges so that a single measure doesn't dominate the search. So a simple least squares distance between the observed data and simulated data may not work well. A better way to handle the situation is by bringing all objectives into a common range (say [0,1]). This can be achieved by dividing each error by the magnitude of its maximum possible error (or just the order of magnitude like 10, 10^2, 10^3...).

Consider two error/distance measures ERR1: Squared distance = ∑(yest-yobs)^2, ERR2: Fraction change = ∑ abs((yest-yobs)/yobs). In ERR2 the magnitude of numerator is proportional to the ERR1, so the difference between them is the denominator. The inverse "yobs" can be thought as weights. If all "yobs" are in similar range then all error terms are weighted similarly. But if there are some outliers with large "yobs", then the function penalizes them with a lower weight. To demonstrate this I will work out a hypothetical example with these two objective functions.

In the figure below, the red dots are the observations noted from a hypothetical experiment which is capturing the relation between two variables X & Y. Theoretically the relation between them is given by Y=X. The experimental observations are vulnerable to external noises and hence the relation is not easily seen. But a careful inspection highlights the outliers at x=3 and x=9. In the code you can see that I added a random normal noise at all points and and extra noise to create outliers at x=3,6,9.

A difference of squared error is sensitive to outliers in the data, designing an objective function incorporating the peculiarities of the case will be more effective 

Now looking at the average error between the family of lines Y=X+intercept and the observations, the lowest minimum error with ERR2 is at intercept=0 whereas the lowest error with ERR1 is at intercept=2. The outlier at x=3 is pulling up the overall fit, proving that the ERR1 is more sensitive to this outlier compared to ERR2.

In the above example the dominant outlier has higher magnitude compared to other observations, hence the weighting 1/yobs reduces the importance attached to these points. So the same weighting scheme may not work in the case where outlier has lower magnitude compared to the general trend. In such a case 1/yobs will add more importance to outliers while fitting. Hence a different approach like multiplication instead of division with yobs will prove useful.

Conclusion :

In this text I discussed the model calibration practiced in a Mathematical model development. To be effective in this phase requires the modeler to assess several items: Mathematical description of relations in the model, A good search strategy and description of error or fitness measure. Here I presented my case on the importance of the design of objective functions with examples to illustrate the point. 

In reality it may be wise to solve the problem on multiple fronts. In my own task I encountered the case of a optimizing the model for multiple objectives simultaneously. I chose to reduce the multiple objectives to a scalar where using fraction error measure helped me avoid worrying over the weights to use. Then I implemented the global optimization strategy Genetic Algorithm. But for the implementation to be practical I vectorized the operations in the program. Further, trying out a different choices of ode solvers gave me an economical option by bringing down the optimization task from an estimated run time of over 20hr+ to slightly over 1hr.

The R code for the last example is hosted on Github

Introduction

In this post I will use Bayesian approach to analysing a dataset containing responses of students to a test.The full code and the dataset can be found on Github link. Bayesian statistics is built on a philosophy different from the classical/frequentist statistics. As a consequence, instead of observing just 1 outcome, we look at multiple possibilities and identify the probability of occurrence of each of these possibilities. The strength of bayesian analysis is in the basic questions we can ask and the simplicity with which the answers can be communicated. We will see such an example in the following

The underlying assumption in the classical way is that "the phenomenon under testing is fixed and the uncertainty observed is because of the experimental and measurement errors". In this approch,

• We start with defining a Null hypothesis and an Alternate hypothesis
• Then we choose an appropriate test to find the p-value
• Based on the decided significance level, α, we can either reject null or fail to reject null

These last 2 statements are usually read as Null hypothesis is wrong and Alternate hypothesis is correct and Null hypothesis is correct and so Alternate hypothesis is wrong respectively. The real world instances may not be so black and white because of the huge number of factors affecting the case (refer chaotic systems).

In Bayesian approach, the underlying phenomenon is believed to be uncertain and this uncertainty reduces as we observe more information. Steps taken to perform a bayesian analysis are

• Define a prior belief for the parameters under study
• Find the posterior distribution of the parameters based on the observed evidence

Posterior is estimated following the bayes rule defined as 

     P(Parameters|Evidence) ∝ P(Evidence|Parameters)*P(Parameters) 

     P(parameters) is the prior distribution 

     P(parameters|Evidence) is the posterior distribution of the parameters conditioned on observing the evidence 

     P(Evidence|parameters) is the likelihood of observing the data given a prior distribution 

The prior distribution can be informative taken from a previous analysis or a good understanding of the system under study. In case where there is no clear understanding of the system, the prior can be non-informative like a uniform distribution. 

Discription of the dataset

In this example we will consider a dataset containing responses of students to a test with Multiple Choice Questions(with a single answer). A subset of the dataset is shown here. Each row has the responses given by a candidate for the corresponding question on the columns. The last column "PointsReceived" is candidate final scoring as calculated by the online judge.

The correct answers to each question and the points alloted to it is shown below.

First, let's quickly check if the scoring calculated by the judge is correct

Looks like the scores calculated are correct hence we can move on to our analysis.

All the questions in this test have single deterministic answer. Hence a candidate response to a question can be either correct or incorrect and so is a bernoulli random variable. Now assuming that all candidates have equal capability to answer all questions, each row of responses are samples from a binomial experiment.

Based on the alloted points, each question can be catergorized as C1(1 pt), C2(2 pts) ans C3(3 pts). This categorization is made by the questioner by choosing unequal weightage to the questions. 

Defining Priors

Lets say a 1pt question can be solved correctly by 80% of the candidates. Since all the candidates are assumed to be equally capable, this means that when posed with this question, the aspirant answers it correctly 80% of the times. Similarly a 2pt question can be expected to be solved correctly by 60% and a 3pt question by 40% of the applicants. Mathematically,

• E(solving correctly | 1pt question) = 0.8
• E(solving correctly | 2pt question) = 0.6
• E(solving correctly | 3pt question) = 0.4

A point to note from above is that E(solving correctly | 1pt question) = 0.8 and not P(solving correctly | 1pt question) = 0.8. This means that there are possibilities(in the multiverse if you may) where a 1pt question can be solved correctly with a probability more or less than 0.8. This is a distinct feature of Bayesian ways as compared to the classical methods.
The distribution of the random variable P(solving a 1pt question correctly) and other two can be taken as beta-distribution. A beta distributed variable has the range [0,1] which is our case. Also, beta-distribution is a conjugate prior for a binomial variable. (ref wikipedia).

There are 2 hyperparameters needed to define the shape of a beta-distribution namely α and β

• mean = α/(α + β)
• concentration = α + β

Mean is the mean value and concentration is to control the uncertainty in the parameter.

Analysis and Inferences

By Fixing the concentration to 50 and the mean values defined before we can calculate the shape parameters. From the observed evidence the posterior hyperparameter values can be found. The prior and posterior parameter values are summarised below (in the format α,β) 

Based on these parameters, samples can be taken from the beta-distribution. The probability density function for probability of solving a question correctly from prior and posterior are shown below.

An interesting observation from the above plots is that, the spread or dispersion in posterior is much less compared to the initial belief. This supports the initial assumption that all the applicants are equally capable to start with and so the uncertainty in the probability of answering correctly is little.

If we now say that any question which is solved by 80% or more of candidates is easy, 40-80% is mediumly tough and <40% is hard question. Then by taking samples from the distribution, we can easily visualise the categorization of a question through a bar graph as shown below.  

Conclusion

In this post I have taken Bayesian approach to analysing a dataset containing responses of applicants to a test. The advantage of using bayesian analysis is in getting an idea of alternate scenarios and their associated chances of occuring. The posterior estimates obtained can be used a prior for the next dataset of similarly capable candidates tested on same questions.

References

While the interest in running the analysis is personal, I had to learn the mathematical concepts somewhere and following references have motivated, inspired and taught me 

• the theory that would not die
• Obscureanalytics blog