Davor Josipovic Just another WordPress blog – rather tryout


Algorithms for detection of drifts and events in air and hydrostatic pressure data

Filed under: Algorithms,Modeling,R,Statistics — Tags: — Davor @ 21:30

Air pressure sensors can start drifting after years of exposure to extreme temperature and weather conditions, producing inaccurate results. The drift is very difficult to detect visually, but relatively easy to detect algorithmically.

This image shows the sensor drifting and the algorithm pointing the most probable start of the drift somewhere at the end of 2013. More information about the algorithm can be found here.

Hydrostatic pressure meter results are susceptible to all kinds of events like systematic tides, temporary effects due to heavy rain, permanent shifts due to equipment or environmental adjustments, and single measurements errors due to sensor inaccuracies. Detecting such events in time series is often tedious and time-consuming.

This image shows how the algorithm decomposes the timeseries into multiple level shifts, two outliers and no temporal changes. More information about the algorithm can be found here.


Calvin’s optimal way home

Filed under: Algorithms,Modeling,Statistics — Tags: , — Davor @ 18:28

This puzzle has been posted here previously but with no correct answer, except for (a). It has also been used by Optiver for the assessment of their new quantitative researchers.

Calvin has to cross several signals when he walks from his home to school. Each of these signals operate independently. They alternate every 80 seconds between green light and red light. At each signal, there is a counter display that tells him how long it will be before the current signal light changes. Calvin has a magic wand which lets him turn a signal from red to green instantaneously. However, this wand comes with limited battery life, so he can use it only for a specified number of times.

a. If the total number of signals is 2 and Calvin can use his magic wand only once, then what is the expected waiting time at the signals when Calvin optimally walks from his home to school?
b. What if the number of signals is 3 and Calvin can use his magic wand only once?
c. Can you write a code that takes as inputs the number of signals and the number of times Calvin can use his magic wand, and outputs the expected waiting time?


My solution to all three questions can be found here: Davor, puzzle solution, 2017. I think it can be interesting to people who have never done statistical modeling and would like to know how it is done.

a. Expected trip time is 8.75 sec. Optimally, wand should be used on red light if the counter is above 20 sec.
b. Expected trip time is 21.32 sec. Optimal wand usage at first light is if the counter is above 31.25 sec, and 20 at the second if there is a charge left.
c. See the solution for the recursive 10-line code that solves the general case. For example, if Calvin has 2 charges and there are 4 lights, then the expected trip time is 11.8 sec with the optimal wand usage.


Extending Berman’s ICE algorithm with spatial information

Filed under: Algorithms,Modeling,Statistics — Tags: , — Davor @ 16:08
Note: this article is a short summary of a larger work I have done here.

Short introduction

Hyperspectral images are like ordinary images, except that they have lots of bands extended beyond the visible spectrum. This extra information is exploited for material identification. For example, in this hyperspectral image a sub-scene is selected – called the Alunite Hill.

The original Cuprite scene from which a 16 x 28 pixel subimage is extracted for analysis.

Subsequently, the ICE algorithm is run which results in the following three material abundance maps contained in matrix \mathbf{W} and material signatures \mathbf{E} called endmembers.

Output of the ICE algorithm: abundance maps for alunite, muscovite and kaolinite minerals, and their hyperspectral signatures, i.e. endmembers.

ICE algorithm

Simplified, the ICE algorithm can be written as an objective function L_{ICE} which measures the error between the actual pixels and the predicted pixels ||\mathbf{X} - \mathbf{W}\mathbf{E}^T ||_F^2 together with a regularization term  V(\mathbf{E}) that measures the size of the simplex formed by the endmembers.

L_{ICE}(\mathbf{E}, \mathbf{W}) = ||\mathbf{X} - \mathbf{W}\mathbf{E}^T ||_F^2 + V(\mathbf{E})

This objective function L_{ICE} is subsequently minimized to get the estimates of the abundance maps and endmembers:

\underset{\mathbf{E}, \mathbf{W}}{argmin} \; L_{ICE}(\mathbf{E}, \mathbf{W})

Spatial information

The idea that spatial information is important stems from the fact that materials in abundance maps are more likely to reflect certain order.

Two abundance maps of the same material. The pixels are identical, but randomly ordered in the right one.

Thus the right abundance map in which the material seems randomly scattered, should be penalized more than the left one where there seems to be a certain order. One way to achieve this is by looking at the adjacent pixels of the abundance maps and see how similar they are. In this specific approach we are calculating the variance of a pixel and its adjacent four pixels. These variances are summed over all pixels of the abundance maps:

S(\mathbf{W})  = \sum_n \sum_m \text{Var}(\mathbf{u}_{n,m})

where \mathbf{u}_{n,m} signifies the vector of the abundance and the adjacent abundances of the n-th pixel and m-th endmember. The new simplified objective functions becomes

L_{ICE-S}(\mathbf{E}, \mathbf{W}) = ||\mathbf{X} - \mathbf{W}\mathbf{E}^T ||_F^2 + V(\mathbf{E}) + S(\mathbf{W})

This new objective function is then minimized:

\underset{\mathbf{E}, \mathbf{W}}{argmin} \; L_{ICE-S}(\mathbf{E}, \mathbf{W})

The abundance maps resulting from the minimization of L_{ICE-S} tend to be more smooth.

Output of ICE-S algorithm. The abundance maps are slightly more smooth when compared with the ICE output above. The central pixel in the third abundance map seems to be completely removed.


My main aspiration here was to give a very concise and simplified version of how the ICE algorithm can be extended with spatial information. The actual topic is much more complex. For those interested: The theoretical foundations, calculus and implementation details can be found here.



Primal Active-Set method for constrained convex quadratic problem

Filed under: Algorithms — Tags: , — Davor @ 10:05

I made an implementation of the Active-Set Method for Convex QP as described in Nocedal, J. e.a., Numerical Optimization, 2ed, 2006, p.472. The code in R can be found on Github. Output with two examples from the book can be found here.

There is an other free package named quadprog that does the same but with the limitation of only accepting positive definite matrices. This can be tweaked to work with positive semi-definite matrices. For example one can convert a positive semi-definite matrix to its nearest definite one with Matrix::nearPD() function. To cope with semi-definiteness in the Active-Set method, Moore-Penrose generalized inverse is used for solving the KKT problem.

Note that the Active-Set method must start with an initial feasible point. Understanding the problem is usually enough to calculate one. Nocedal describes a generic “Phase I” approach for finding a feasible point on p.473.

Powered by WordPress