Showing posts with label statistics. Show all posts
Showing posts with label statistics. Show all posts

Saturday, September 20, 2025

Book Review: Statistics every Programmer Needs

I recently read Statistics every Programmer Needs by Gary Sutton. I am probably a good target audience for the book since I used to be a software developer that transitioned into data science some 10 years ago, then into machine learning with neural networks and transformers, and more recently, to Generative AI with Large Language Models. During this time, I have read numerous books on statistics in an effort to pick up what I didn’t know (being largely self-taught, there is plenty I didn't and still don't know). I think this book stands out not only as a thorough and practical introduction to statistics, but also provides coverage to areas one would normally consider peripheral to statistics but still useful in practical data science scenarios, such as Linear Programming, PERT/CPM, etc.

The book takes a very hands-on approach to each area, starting with business problems often faced by programmers, and outlines how statistical techniques (pertinent to that area) can be used to address these problems. It starts with foundational concepts but goes on to cover advanced concepts across statistics, machine learning, optimization, and project management. The book is organized into the following 14 chapters.

The Foundation (Chapter 1) begins by laying a solid groundwork. Readers are introduced to core statistical concepts, both descriptive (mean, mode, median) and inferential (confidence intervals, p-values), ensuring they grasp the basics before progressing. The inclusion of regression, optimization, simulation, and machine learning in the foundational chapter sets the tone for the book’s broad scope.

Probability and Counting Principles (Chapter 2) covers continuous and discrete variables and how they differ, permutations, combinations, and key probability functions (PDF, PMF, CDF). What was interesting for me is how permutations and combinations are described using basic probability concepts.

Probability Distributions (Chapter 3) covers the essential probability distributions—Gaussian, Binomial, Uniform, and Poisson. This chapter also covers conditional probability and Bayes’ rule with various practical applications.

Chapters 4 & 5 cover Linear and Logistic Regression respectively, Bonus material here (which I didn’t expect to see at least) were discussions around data normalization, residual analysis and multi-collinearity. Model evaluation is covered in depth for both varieties of models, as well as discussion of popular metrics used to evaluate these models.

Chapter 6 covers Decision Trees and Random Forests, the next major category of traditional ML models. The book has a solid introduction to decision trees and random forests, including how to interpret feature importance and use GINI impurity measures. I had hoped for some coverage of Gradient Boosted Trees since we were already discussing trees, but maybe that will come in the next edition.

Time Series Analysis (Chapter 7) is tackled with impressive depth, usually I would expect this subject to have its own book. However, the author does a good job of providing a good useful introduction to Time Series – covering forecasting, ARIMA models, exponential smoothing, stationarity testing (including the Augmented Dickey-Fuller test), trends, and seasonality. The chapter’s coverage of ACF/PACF plots and different exponential smoothing models (SES, DES, Holt-Winters) is thorough, making it a valuable reference for people working with temporal data and autoregressive models.

Chapter 8 covers Optimization using Linear Programming, an area I would expect to see covered in a book on Operations Research rather than Statistics. But the coverage is practical and complete, focusing on modeling business problems as optimization problems and solving them using Linear Programming libraries provided by scipy.optimize.

Chapter 9 covers Simulation using Monte Carlo techniques. As before, not something I would have expected in a Statistics book, but definitely a useful tool to have in one’s Data Science toolbox. As with the other chapters, multiple business scenarios are described and modeled with probability distributions, and Monte Carlo simulations performed on them to elicit useful insights.

Decision Methods and Markov Analysis (Chapters 10 & 11) cover Decision-making frameworks (maximax, maximin, minimax regret, expected value decision trees) and Markov analysis (transition probabilities, equilibrium, and absorbing states). Taken together, they could serve as a gateway for deeper explorations into Bayesian Networks and other Probabilistic Graphical Models.

The chapter on Benford’s Law (Chapter 12) for fraud detection is another unique touch, introducing readers to mantissa statistics. So is the chapter on Project Management (Chapter 13), which presents quantitative methods in project management (WBS, PERT, CPM, critical path)with actionable insights, bridging the gap between theory and project execution.

The concluding chapter on Statistical Quality Control (Chapter 14) is packed with practical content—control charts (p, np, c, g, etc.), UCL/LCL, and key metrics—making it invaluable for readers in manufacturing, operations, or quality assurance roles.

I thought that the book is ambitious in scope but succeeds in providing both breadth and depth, managing to hit all the high points without impacting the quality of each. As I mentioned earlier, its coverage goes beyond just statistics, making it a bargain since you get to learn useful statistics and quantitative techniques from a single book. I found both areas to be described in a very hands-on, example driven manner,often highlighting concepts and metrics that are overlooked in more traditional texts, thus making it a useful reference for software professionals (DS and non-DS alike).

Saturday, February 03, 2024

Book Report: Allen B Downey's Probably Overthinking It

I have read Allen Downey's books on statistics in the past, when trying to turn myself from a Software Engineer into what Josh Wills says a Data Scientist is -- someone who is better at statistics than a Software Engineer and better at software than a statistician (with somewhat limited success in the first area, I will hasten to add). Last year, I had the good fortune to present at PyData Global 2023 (the video is out finally!) so had a free ticket to attend, and one of the talks I really enjoyed there was Allen Downey's talk Extremes, Outliers and GOATs: on life in a lognormal world. In it, he mentions that this is essentially the material from Chapter 4 of his book Probably Overthinking It. I liked his talk enough to buy the book, and I wanted to share my understanding of this book with you all, hence this post.

The book is not as dense as a "real" book on stats like say The Elements of Statistical Learning but is definitely not light reading. I tried reading it on a flight from San Francisco to Philadelphia (and back) and found it pretty heavy going. While the writing is lucid and illustrated with tons of well-explained and easy to understand examples, most of these were new concepts to me, and I wished I took notes after each chapter so I could relate all these concepts together enough to reason about them rather than just learn about them. So I did another pass through the book, this time with pen and paper, and I now feel more confident about talking to other people about it. Hopefully, this is also helpful for folks who have done (or planning to do) the first pass on the book but not the second.

Most people who are new to statistics (me included) lay great store in the Gaussian (Normal) distribution to explain or model various datasets. Chapter 1 challenges this idea and demonstrate that while individual traits may follow a Gaussian distribution, a combination of such traits can be a very restrictive filter. In other words, almost all of us are weird (i.e. not normal). For me, it also introduces the Cumulative Distribution Function (CDF) as a modeling tool.

The second chapter introduces the Inspection Paradox, which explains why it always seems like our wait time for the next train is longer then the average wait time between trains, among other things. The explanation lies in the sampling strategy -- if we sample our data from the population, we may get a skew from oversampling from over-represented populations. It also describes a practical use case of this paradox to detect COVID superspreaders.

The third chapter describes what the author calls Preston's paradox, based on a 1976 paper by Samuel Preston. The paradox is that even if every woman has fewer children than her mother, the average family size can increase over time. The paradox is explained by an idea similar to the Inspection Paradox, i.e. because there are more women in existence from large families than small ones, a larger proportion of women would end up having large families than small ones, and overall that contributes to an increase in family size. The opposite can hold true as well, as demonstrateed by the loosening of reproductive restrictions in China in the aftermath of China's one-child policy not having the desired effect in boosting family sizes.

Chapter 4 is the one the author talked about in the PyData Global talk. In it, he demonstrates that certain attributes are better explained by a log-normal distribution, i.e. taking the log of the values in the distribution, rather than our familiar Gaussian distribution. This is especially true for outlier type distributions, such as performance numbers of GOAT (Greatest Of All Time) athletes compared to the general population. The explanation for this is that GOAT performance is almost always a multiplicative combination of innate human prowess (nature) and these skills being effectively harnessed and trained (nurture) plus a whole lot of other factors that all have to line up just so for the event to happen, and whose contributions to the target are therefore multiplicative rather than additive, hence the effectiveness of the log-normal distribution over the normal one.

Chapter 5 explores different survival characterstics of different populations and classifies them as either NBUE (New Better than Used in Expectation) and NWUE (New Worse than Used in Expectation). The former would apply for predicting the remaining life of lightbulbs with use, and the latter would apply for predicting cancer survivability and child mortality over time. Using child mortality statistics, the author shows that as healthcare improves and becomes more predictable across age categories, the NWUE distribution changes to resemble more closely a NBUE distribution.

Chapter 6 explores Berkson's Paradox, where a sub-sample selected from a population using some selection criteria can create correlations that did not exist in the population, or correlations that are opposite to that observed in the population. Berkson originally pointed out the paradox as a warning about using hospital data (sub-sample) to make conclusions about the general population. The selection criteria restrict the general population in specific ways, leading to a change in composition of the traits in the sub-sample, thus leading to the paradox.

Chapter 7 warns about the dangers of interpreting correlation as causation, something most of us have probably read or heard about many many times in the popular Data Science literature. The main case study here are moms who smoke (or don't smoke) and their low birth weight (LBW) babies. A study concluded that while smoker's were more likely to give birth to LBW babies, and LBW babies had a higher mortality rate, the mortality rate of LBW babies whose mothers smoked was 48% lower than those whose mothers didn't smoke. Further LBW babies of non-smokers also had higher rate of birth defects. Interpreting this correlation as causation, i.e. not heeding the warning, it seems like maternal smoking is beneficial for LBW babies, protecting them from mortality and birth defects. The explanation is that maternal smoking is not the only cause of LBW babies, and birth defects may be congenital and not linked to smoking. These two factors mean that there are biological explanations for LBW other than maternal smoking. This and a few other examples segue naturally into a brief high-level introduction to Causal Reasoning, which I also found useful.

Following on from GOAT events being better represented by log-normal rather than normal distributions, Chapter 8 describes applying this to model extremely rare events (such as earthquakes and stock market crashes), and concludes that while the log-normal distribution is more "long-tailed" than a Gaussian, rare events have an even longer tail that is better modeled by log-Student-t (or Log-t) distibution (Student-t is a Gaussian with longer / fatter tails). It also introduces the idea of a Tail distribution (the inverse of a CDF, a survival chart is a tail distribution chart). The author also makes a brief reference to Nassim Taleb's Black Swan events, saying that the ability to model and predict them make them more of Gray Swans.

Chapter 9 talks about the challenges in ensuring algorithmic fairness to all recipients of its predictions, which is very relevant given the many paradoxes the book has already covered. In this chapter, the author describes Bayes rule without mentioning it by name, calling it the "base rate" and the difference between the prior and posterior probabilities the "base rate fallacy". He also covers other aspects of fairness, citing differences across groups that an algorithm often does not see. This last part seemed to me to be related to the Inspection Paradox described earlier in the book.

Chapter 10 describes Simpson's Paradox, where sub-populations can exhibit similar correlations across the sub-populations but where the same traits are anti-correlated in the conbined population. To some extent, this seems related to Berkson's law. Among the examples cited, there is one about penguins, where within each species, the beak size and body size are correlated, but across species, they are anti-correlated. The explanation here is that there is a biological reason for the correlation within the species, but the anti-correlation is just a statistical artifact (correlation != causation in action I guess?).

Chapter 11 is about how certain instances of Simpson's Paradox can be explained as a combination of other underlying factors. It is a trusim that people get more conservative as they get older (i.e. if you are not a liberal when you are young, you have no heart, and if you are not a conservative when old, you have no brain). However, within each age group, it is observed that people actually get more liberal over time. This is explained as a combination of the age effect, the period effect, and the cohort effect. The age effect shows a positive correlation between adherence to traditional beliefs (conservativeness) and age. However, within each age group, it is observed that people get more liberal over time, i.e. the cohort effect. Finally the period effect deals with specific events during the time period under consideration, and this covers older people dying out and being replaced with younger (and more liberal) people.

Chaoter 12 continues the discussion from the previous chapter and brings in the idea of the Overton Window, which dictates what views are considered acceptable at any particular point in time, and which changes over time as well. So what was thought to be liberal in decades past is now considered more conservative. So while an individual may get more liberal with time, the Overtom Window has shifted faster towards liberalism. This can explain why an individual may find themselves getting more conservative as they age, relative to the world around them.

Overall, I enjoyed this book. I think the most impressive thing about this book was its use of generally available datasets to model physical and social environments, and using simulations to control for certain aspects of these data experiments. Also, I think I learned a few things about corner cases in Statistics which I think may be useful when reasoning about them in future. I hope I have sparked your curiosity about this book as well.

Saturday, June 06, 2020

Bayesian Analysis with Python Exercise Solutions


In my quest to learn more about Bayesian Analysis, I spent the last few weeks working through the book Bayesian Analysis with Python - 2nd Edition by Osvaldo Martin. In their great wisdom, Amazon now only allows reviews by people who have been spending "real" money on Amazon, as opposed to those like me using "fake" money in the form of gift cards and credit card rewards. So I guess I will have to hold off posting my review of this great book until my fake money runs out. Meanwhile, you can read other people's review on Amazon. In this post, I would like to share my notebooks where I have worked out the exercises for Chapters 2-7 (there are 9 chapters, the first one is an introductory one, and the last two chapters cover areas only indirectly related to Bayesian Analysis, or more specifically, Bayesian analysis using Monte Carlo Markov Chain (MCMC) techniques. Here are my notebooks on Github. No guarantees, of course, since I am learning as well.


The author Osvaldo Martin is one of the developers of the ArviZ project, an exploratory analysis tool for Bayesian models. The ArviZ project was spun-off from the PyMC3 project, and many PyMC3 calls such as pm.traceplot() are actually calls to az.plot_trace() under the hood. According to the ArviZ website, it also supplies functionality for other Bayesian libraries, such as PyStan, Pyro, and TF Probability. Martin is also a PyMC3 maintainer, and the book reflects that. Compared with my previous foray into PyMC3, where I was trying to convert JAGS models taught in the two Coursera courses from the University of California, Santa Cruz (UCSC), into equivalent PyMC3 models using online documentations and other Internet resources, this time the focus was on exploring various code features offered by PyMC3 and ArviZ.

The book is billed as an introductory one on the PyMC3 site, but I thought it had good followup material for someone who has already taken the two UCSC Coursera courses. The focus on those courses was on the theory, to explain why you would choose one distribution over the other, linear and logistic regression, and some exposure to hierarchical and mixture models. This book contains a more code-heavy exposition of all these subjects, plus an introduction to Dirichlet processes as part of the Mixture models chapter, and an introduction to Gaussian processes.

Of course, I still have much more to learn about all this. Just today, I came across a link to Nicole Carlson's PyMC3 talk at PyData Chicago 2016, where she walks through some toy problems with PyMC2. She introduces the ideas of using Theano shared variables for input and output here, so you can use the model in a more Scikit-Learn like way, training with one dataset and evaluating with another. She also covers the newer PyMC3 Variational Inference API (ADVI), an alternative to its older MCMC API, and how to use the output of the ADVI model as input to an MCMC sampling model. In addition, she also describes serializing a PyMC3 model by pickling the trace and using it as a generative model. I also found the notebooks for her subsequent talk at PyData Chicago 2017, where she describes her work to wrap PyMC3 models in a more familiar Scikit-Learn like interface.

Another recommendation from the PyMC3 site, as well as from this book, was the "puppy book" Doing Bayesian Data Analysis by John Kruschke. Fortunately, that is available in my employer's online library, so I guess that is another book I won't buy from Amazon with fake or real money.

Thats pretty much all I had for today. Obviously the best way to learn this stuff is to try doing it yourself, but sharing these notebooks on Github in case someone gets stuck and can benefit from them. In case you do end up using them, please let me know what you think.

Sunday, April 05, 2020

Experiments with COVID-19 Patient Data


An interesting (and beneficial, in spite of the otherwise less than ideal situation) side-effect of the COVID-19 pandemic has been that many organizations, both commercial and academic, are coming together looking for ways in which they can work together to eradicate the disease. Many of the collaborations involve sharing datasets, the most famous of which is the COVID-19 Open Research Dataset (CORD-19), a collection of 47,000 (and growing) scientific papers about COVID-19. Some others are offering the CORD-19 dataset processed through their pipelines, or hosted using their products (for example, graph database or search engine). Some are holding seminars, and sharing their expertise with the rest of the world. At Elsevier, we have a grassroots Data Science team of more than 225 employees, looking at the problem from different disciplines and angles, and working to find solutions to address the crisis. The LinkedIn article Elsevier models for COVID19 bio-molecular mechanisms describes some contributions that were driven by work from this team using one of our tools, and hopefully there will be more soon. In addition, about 46% of the papers in the CORD-19 dataset come from Elsevier, and we are looking at ways of making more available.

In the spirit of learning everything I could about COVID-19, I attended the day-long COVID-19 and AI: A Virtual Conference organized by the Stanford Human-AI (HAI) group. One of the speakers was Prof. Nigam Shah, who spoke about his Medical Center's Data Science Response to the Pandemic, and described the types of Data Science models that can inform policy to combat the virus. In addition, he also wrote this Medium post about Profiling presenting symptoms of patients screened for SARS-Cov-2 where he used the same diagram for his unified model, which is what caught my eye. Hat tip to my colleague Helena Deus for finding and posting the link to the article on our internal Slack channel.

In any case, the Medium post describes a text processing pipeline designed by Prof. Nigam's group to extract clinical observations from notes written by care providers at the Emergency Department of Stanford Health Care, when screening patients for COVID-19. The pipeline is built using what look like rules based on the NegEx algorithm among other things, and Snorkel to train models that recognize these observations in text using these noisy rules. The frequency of these observations were then tabulated and probabilities calculated, ultimately leading to an Excel spreadsheet, which Prof. Nigam and his team were kind enough to share with the world.

There were 895 patients considered for this dataset, of which 64 tested positive for SARS-Cov-2 (new name is COVID-19) and 831 tested negative. So at this point in time, the prevalence of COVID-19 in the cohort (and by extension, possibly in the broader community) was 7.2%. The observations considered in the model were the ones that occurred at least 10 times across all the patient notes.

So what can we do with this data? My first thought was a symptom checker, which would compute the probability that a particular patient test positive given one or more of the observations (or symptoms, although I am using the term a bit loosely, there are quite a few observations here that are not symptoms). For example, if we wanted to compute the probability of the patient testing positive given that the patient exhibits only cough and no other symptom, we would denote this as P(D=True|S0=True, S1=False, ..., S49=False).

Of course, this depends on the simplifying (and very likely incorrect) assumption that the observations are independent, i.e., the fact that a patient has a cough is independent from the fact that he has a sore throat. Also, the other thing to remember is that predictions from the symptom checker will be dependent on the correct value of the current disease prevalence rate. The 7.2% value we have is only correct for the time and place where the data was collected, so will need to be updated accordingly if we wish to use the checker even with all its limitations. Here is a schematic of the model.


Implementation wise, I initially considered a Bayesian Network, using SQL tables to model it as taught by Prof. Gautam Shroff in his now-defunct Web Intelligence and Big Data course on Coursera (here's a quick note on how to use SQL tables to model Bayesian Networks since the technique, even though its super cool, does not appear to be mainstream), but I realized (thanks to this Math StackExchange discussion on expressing Conditional Probability given multiple independent events), that the formulation can be much more straightforward, as shown below, so I used this instead.


The idea of using the proportionality relationship is to normalize the numerator by computing P(D=True|∩Sk).P(D=True) and P(D=False|∩Sk).P(D=False) and divide by the sum to get the probability of a positive test given a set of symptoms. Once that was done, it led to several additional interesting questions. First, what happens to the probability as we add more and more symptoms? Second, what happens to the probability with different prevalence rates? Finally, what is the "symptom profile" for a typical COVID-19 patient based on the data? Answers to these questions and the code to get to these answers can be found in my Github Gist here.

I have said it before, and given that people might tend to grasp at straws because of the pandemic situation, I am going to say it again. This is just a model, and very likely an imperfect one. Conclusions from such models are not a substitute for medical advice. I know most of you realize this already, but just in case, please do not use the conclusions from this model to make any real life decisions without independent verification.

Saturday, March 28, 2020

Python notebooks for Bayesian Analysis Courses on Coursera


I recently completed the Coursera courses Bayesian Statistics: From Concept to Data Analysis and Bauesian Statistics: Techniques and Models, taught by Prof. Herbert Lee and Mathew Heiner of the University of California, Santa Cruz. I did both in audit mode, so "completed" is not totally accurate, since the second course did not allow submission of quiz answers without paying for the course. But the content for both are free and excellent, and I learned a lot from them, and highly recommend them if you are interested in the subject of Bayesian Analysis. Please be aware that the courses are somewhat elementary, it was a good way for someone like me, curious about but not very knowledgeable about Bayesian Analysis, to get to a point where I can hopefully explore the subject on my own. So if you are like me, you will find the courses useful, otherwise probably not.

Both courses use R as the programming language. The first course is more math and less programming, but covers concepts which are essential in the second course. In fact, I started on the second course because I was curious about MCMC (Markov Chain Monte Carlo), but found myself out of my depth within the first week, so I ended up having to do Prof. Lee's course first. And even though it's called Concepts to Data Analysis, this is much more than your high school statistics course (my level before I took the course, give or take). It starts with probability concepts, then goes off into different kinds of distributions and when you should use them, how to do inference with these distributions and Bayes theorem, both for discrete and continuous data. At the end of the course, you will know which distributions to use when, and what to look for when trying to draw conclusions from a given distribution. It also covers linear regression, both single and multiple, from a statistician's rather than a machine learning perspective.

The second course is taught by Mathew Heiner, a doctoral student at UCSC. This expands on Prof. Lee's course, starting with simple conjugate models (this is where I realized I was out of my depth the first time around, BTW) and moving on to MCMC models for binary, continuous, and count data, as well as how to compose them into hierarchical models to account for uncertainty in our knowledge. It also covers Metropolis-Hastings and Gibbs sampling methods. The course is very example-driven, using small datasets included in the R platform, to explain each concept. The MCMC library used in the course is rjags, which depends on JAGS (Just Another Gibbs Sampler).

Going into the course, I had some understanding (superficial in hindsight) of MCMC, and my main motivation was to learn enough theory to work intelligently with PyMC3, a Python toolkit for Probabilistic Programming. I figured that going through the course to learn the theory and reimplementing the R and JAGS examples in Python and PyMC3 will allow me to learn both faster (kind of how joint learning sometimes works better in Machine Learning models), so that's what I did. These are modeled as Jupyter notebooks with short text annotations, code, and outputs, so you can read them if you like, but you will probably benefit more from doing the exercises yourself and using my notebooks as cheat sheets for when you are stuck. All notebooks are runnable without any additional data. The examples used datasets built into the R platform, which Vincent Arel-Bundock has been kind enough to package and host at his R-Datasets repository. My notebooks automatically pull the data from his repository if they are not already downloaded.

The notebooks are in my sujitpal/bayesian-stats-examples Github repository, each course is in its own subfolder. Direct links to notebooks for each course are provided below for convenience as well, hopefully you find them useful as you navigate your way around the world of Bayesian analysis and PyMC3.


I am currently exploring this subject a bit more using the book Bayesian Analysis with Python 2/ed by Osvaldo Martin. The book is recommended from the PyMC3 github page, and so far, I find it covers the Coursera course material, and then some, even though it is listed as an introductory book. The PyMC3 Tutorial is also an excellent resource, and I have used it as a reference when reimplementing JAGS models from the course. The book also mentions the Arviz package for exploratory analysis of Bayesian models, which is part of the effort around the move to PyMC4 (see below), and is being led by the author.

Another book I want to mention is the one where I first learned about PyMC3 -- Bayesian Methods for Hackers by Cam Davidson-Pilon. The book is fantastic and only around 250 pages, and contains many code examples and graphs. It resembles a series of very well-written Jupyter notebooks, which is how Davidson-Pilon has effectively also provided an open source version of the book on Github. But I found it very dense the first time I read it, probably because I didn't have sufficient background in statistics to follow it through its mental leaps. In a subsequent partial re-read after completing the courses above, I did end up with an easier read.

Finally, I wanted to address a concern that I had, and I think many others might have too. PyMC3 depends on Theano for its fast numerical computing backend, and the first time I looked it, I learned that LISA Labs, the group at the University of Montreal that created and maintained Theano, had decided that it was time to move on from Theano and discontinued support. At around the same time, the PyMC4 project was born, and its objective was to provide a PyMC3 like API on top of the Tensorflow Probability library. At the time, the future of PyMC3 seemed uncertain, and I figured it might be safer to wait until PyMC4 became available, rather than spending time learning PyMC3 and having my newly acquired skills go obsolete soon after. However, 6-10 months since then, PyMC4 is still pre-release, and the PyMC3 team has committed to supporting Theano as it relates to PyMC3, so I have more confidence that the effort to learn PyMC3 will not go to waste. The article Theano, Tensorflow, and the future of PyMC3 posted by Chris Fonnesbeck, creator of PyMC3, provides more detail around this.

I hope the notebooks are useful. In keeping with Coursera terms of use, I have not published notebooks containing quiz answers, even though their usefulness solely for quiz answers is doubtful. This is because because the Python/PyMC3 models sometimes produce slightly different results from their R/JAGS counterparts described in the course, probably because of numerical precision and algorithm differences. On the other hand, the models on which the quiz questions are based are sometimes interesting because they illustrate concepts mentioned in the classes, so being able to publish them would probably have been helpful.


Monday, February 25, 2019

Autoencoder as a Feature Detector for Complex Classifiers


Using Autoencoders as anomaly detectors is not a new idea. For examples, see here, here and here (results of a google search "autoencoder anomaly detection"). However, one other possibility could be to use an autoencoder as a feature detector for specific features, for example, to look for specific kinds of objects in images. So assuming you wanted to detect domestic animals in images, you could train a set of autoencoders, each trained on a single animal, then run the image through the set. The intuition is that, assuming the image has a cat in it, an autoencoder trained on cat images should be able to reconstruct the image with lower reconstruction error than it would for other animals. In this post, I describe an experiment where I train an autoencoder on cat images from the Kaggle Cats and Dogs Dataset, and compare reconstruction errors on a held-out set of dog and cat images.

The Kaggle Cats and Dogs dataset contains images of 12,500 cats and 12,500 dogs, both color and monochrome. I picked a random sample of 2,500 images each of cats and dogs as my held-out test set, and trained my cat autoencoder on the remaining 10,000 cat images. The autoencoder was trained on monochrome images for 50 epochs. I then compared a few random images and the reconstructed image (input and output of the autoencoder) and it looks like the autoencoder has learned to reconstruct the images quite nicely, as you can see below.






As expected, the autoencoder seems somewhat better at reconstructing cat images compared to dog images (Root Mean Square Error (RMSE) between original and reconstructed image for the cat and the dog image above is 0.04304 and 0.07650 respectively). I then computed reconstruction errors (RMSE again) for the held out set of 2,500 cat and dog images and plotted the distribution of reconstruction errors, which is shown below.


As you can see, while there is quite a bit of overlap between the two distributions, the distribution of reconstruction errors for cat images peaks earlier and is skewed left a little more than the distribution of reconstruction errors for dog images. It also looks like the distributions are approximately normal, so we can use the t-test to detect if there is statistically significant difference between them, as explained on this page. Turns out that the t-value is 16.94 ≫ 3.291, so we can conclude that the distributions are different at 99% confidence level.

Since the distributions of reconstruction errors are approximately normal, and since we can represent a normal distribution by the mean and standard deviation, we can now compute the probability of an image being that of a cat based on its reconstruction error. For example, our cat and dog exampleis have a 80.3% and 0.15% chance of being a cat respectively (code for calculating these probabilities is provided later in the post). This gives us a way to use the autoencoder as a binary cat classifier.

So that was a brief overview of the idea and the results from my experiments. Now I will describe the experiment in more detail, along with the code, so you can replicate the results for yourself, should you find it interesting enough to do so.

Data cleaning


The images are provided in a directory structure that looks like this. The Cat and Dog directories each contain 12,500 images of varying sizes, both monochrome and color. We randomly sample 2500 images each from each collection as our held-out test set and write the filenames to flat files test-cats.txt and test-dogs.txt. The remaining 10,000 cat images are written into train-cats.txt. We use these files as inputs to a custom image generator.

1
2
3
4
5
6
7
8
PetImages/
|
+-- Cat/
|   |
|   +-- cat images
+-- Dog/
|   |
|   +-- dog images

We standardize our input images to all be monochrome and size 128 x 128 pixels. Most images have 3 channels (Red, Green, Blue) which we coalesce to a single channel using the Luminosity formula defined here. Some images had 4 channels (RGB + Alpha channel), for which we first discarded the Alpha channel and applied the same Luminosity formula to the remaining channels. A few files could not be read and were discarded. A few others had issues with its EXIF headers, which seems a spurious problem without a solution in the Pillow package according to this issue, which I ended up discarding as well. Because we are referencing the image files via text files, "discarding" an image is just commenting out its file name. In all, 24 filenames in all were discarded in this manner, 13 from the training set and 5 and 6 from the cat and dog test sets respectively.

The code to pre-process incoming images into a (128, 128) monochrome images is shown below. One thing to note is that the skimage.transform.resize() method not only resizes the image, but also silently rescales the pixel intensities from the [0, 256] range to the [0, 1] range.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
from skimage.transform import resize

import os
import matplotlib.pyplot as plt
import numpy as np

IMAGE_HEIGHT, IMAGE_WIDTH, NUM_CHANNELS = 128, 128, 1

def get_processed_image(image):
    image_r = resize(image, (IMAGE_HEIGHT, IMAGE_WIDTH))
    if len(image_r.shape) == 2:  # 2D image
        image_r = np.expand_dims(image_r, axis=2)
    else:
        if image_r.shape[2] == 3:
            image_r = (image_r[:, :, 0] * 0.21 + 
                       image_r[:, :, 1] * 0.72 + 
                       image_r[:, :, 2] * 0.07)
        elif image_r.shape[2] == 4:
            image_r = image_r[:, :, 0:3]
            image_r = (image_r[:, :, 0] * 0.21 + 
                       image_r[:, :, 1] * 0.72 + 
                       image_r[:, :, 2] * 0.07)
        else:
            pass
        image_r = np.expand_dims(image_r, axis=2)
    assert(image_r.shape == (IMAGE_HEIGHT, IMAGE_WIDTH, NUM_CHANNELS))
    return image_r

Discarding the problem images was semi-manual. I would repeatedly loop through the three lists of files, looking for exceptions, using the following code. Once I found one, I would comment it out so it would get ignored in the next run. To save time, I did a little trick with global variables, so I could start with the file I left off, that bit is commented out. The code below should run without errors once all the offending files are commented out of the list.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
import warnings

warnings.filterwarnings("ignore", ".*will be changed.*")

DATA_DIR = "../data"

CAT_IMAGE_DIR = os.path.join(DATA_DIR, "PetImages", "Cat")
DOG_IMAGE_DIR = os.path.join(DATA_DIR, "PetImages", "Dog")

CAT_TRAIN_FILE = os.path.join(DATA_DIR, "train-cats.txt")
CAT_TEST_FILE = os.path.join(DATA_DIR, "test-cats.txt")
DOG_TEST_FILE = os.path.join(DATA_DIR, "test-dogs.txt")

#LAST_FOUND_JPG = "573.jpg"
#start_processing = False

def get_source_image(imagedir, imagefile):
    image = plt.imread(os.path.join(imagedir, imagefile))
    return image

def check_images(imagedir, imagelist):
#    global start_processing
    num_processed = 0
    legend = os.path.basename(imagelist)
    flist = open(imagelist, "r")
    for line in flist:
        if line.startswith("#"):
#            if line.strip() == "#" + LAST_FOUND_JPG:
#                start_processing = True
            continue
        if num_processed % 100 == 0:
            print("{:s}: {:d} lines loaded...".format(legend, num_processed))
#        if not start_processing:
#            num_processed += 1
#            continue
        imagefile = line.strip()
        try:
            image = get_source_image(imagedir, imagefile)
            image_p = get_processed_image(image)
            if len(image_p.shape) != 3:
                print(imagefile, image.shape, image_p.shape)
        except IndexError as e1:
            print(imagefile, image.shape)
            raise e1 
        except OSError as e2:
            print(imagefile)
            raise e2
        except UserWarning as e3:
            print(imagefile)
            raise e3
        except AssertionError as e4:
            print(imagefile)
            raise e4
        num_processed += 1
    
    flist.close()
    print("{:s}: {:d} lines loaded, COMPLETE".format(legend, num_processed))

check_images(CAT_IMAGE_DIR, CAT_TRAIN_FILE)
check_images(CAT_IMAGE_DIR, CAT_TEST_FILE)
check_images(DOG_IMAGE_DIR, DOG_TEST_FILE)

Model Training


The autoencoder design is based on the one Aditya Sharma describes in his blog post on Reconstructing Brain MRI Images using Deep Learning (Convolutional Autoencoder). The network consists of 2 blocks of convolution and maxpooling layers for the encoder, followed by a final convolutional block, followed by 2 blocks of convolution and upsampling for the decoder. Input to the network is mini-batches (of 128 images each), each image having the shape (128, 128, 1). Here is the code to build the autoencoder. The loss function is Mean Square Error, since we want to minimize pixel-wise reconstruction error, and the optimizer we use is RMSProp with a learning rate of 0.0003.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
import tensorflow as tf

def build_model_2():
    input_img = tf.keras.Input(shape=(128, 128, 1))
    # encoder
    conv_1 = tf.keras.layers.Conv2D(32, (3, 3), activation="relu",
                                    padding="same")(input_img)
    conv_1 = tf.keras.layers.BatchNormalization()(conv_1)
    conv_1 = tf.keras.layers.Conv2D(32, (3, 3), activation="relu",
                                    padding="same")(conv_1)
    conv_1 = tf.keras.layers.BatchNormalization()(conv_1)
    pool_1 = tf.keras.layers.MaxPool2D((2, 2))(conv_1) # (64, 64, 32)
    
    conv_2 = tf.keras.layers.Conv2D(64, (3, 3), activation="relu",
                                    padding="same")(pool_1)
    conv_2 = tf.keras.layers.BatchNormalization()(conv_2)
    conv_2 = tf.keras.layers.Conv2D(64, (3, 3), activation="relu",
                                    padding="same")(conv_2)
    conv_2 = tf.keras.layers.BatchNormalization()(conv_2)
    pool_2 = tf.keras.layers.MaxPool2D((2, 2))(conv_2) # (32, 32, 64)
    
    conv_3 = tf.keras.layers.Conv2D(128, (3, 3), activation="relu",
                                    padding="same")(pool_2)
    conv_3 = tf.keras.layers.BatchNormalization()(conv_3)
    conv_3 = tf.keras.layers.Conv2D(128, (3, 3), activation="relu",
                                    padding="same")(conv_3)
    conv_3 = tf.keras.layers.BatchNormalization()(conv_3) # (32, 32, 128)
    
    # decoder
    conv_4 = tf.keras.layers.Conv2D(64, (3, 3), activation="relu",
                                    padding="same")(conv_3)
    conv_4 = tf.keras.layers.BatchNormalization()(conv_4)
    conv_4 = tf.keras.layers.Conv2D(64, (3, 3), activation="relu",
                                    padding="same")(conv_4)
    conv_4 = tf.keras.layers.BatchNormalization()(conv_4)
    up_1 = tf.keras.layers.UpSampling2D((2, 2))(conv_4) # (64, 64, 64)
    
    conv_5 = tf.keras.layers.Conv2D(32, (3, 3), activation="relu",
                                    padding="same")(up_1)
    conv_5 = tf.keras.layers.BatchNormalization()(conv_5)
    conv_5 = tf.keras.layers.Conv2D(32, (3, 3), activation="relu",
                                    padding="same")(conv_5)
    conv_5 = tf.keras.layers.BatchNormalization()(conv_5)
    up_2 = tf.keras.layers.UpSampling2D((2, 2))(conv_5) # (128, 128, 32)
    
    decoded = tf.keras.layers.Conv2D(1, (3, 3), activation="sigmoid",
                                     padding="same")(up_2) # (128, 128, 1)
    
    autoencoder = tf.keras.Model(inputs=[input_img], outputs=[decoded])
    rmsprop = tf.keras.optimizers.RMSprop(lr=3e-4)
    autoencoder.compile(optimizer=rmsprop, loss="mse")
    return autoencoder

Or in tabular form.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
input_1 (InputLayer)         (None, 128, 128, 1)       0         
_________________________________________________________________
conv2d (Conv2D)              (None, 128, 128, 32)      320       
_________________________________________________________________
batch_normalization (BatchNo (None, 128, 128, 32)      128       
_________________________________________________________________
conv2d_1 (Conv2D)            (None, 128, 128, 32)      9248      
_________________________________________________________________
batch_normalization_1 (Batch (None, 128, 128, 32)      128       
_________________________________________________________________
max_pooling2d (MaxPooling2D) (None, 64, 64, 32)        0         
_________________________________________________________________
conv2d_2 (Conv2D)            (None, 64, 64, 64)        18496     
_________________________________________________________________
batch_normalization_2 (Batch (None, 64, 64, 64)        256       
_________________________________________________________________
conv2d_3 (Conv2D)            (None, 64, 64, 64)        36928     
_________________________________________________________________
batch_normalization_3 (Batch (None, 64, 64, 64)        256       
_________________________________________________________________
max_pooling2d_1 (MaxPooling2 (None, 32, 32, 64)        0         
_________________________________________________________________
conv2d_4 (Conv2D)            (None, 32, 32, 128)       73856     
_________________________________________________________________
batch_normalization_4 (Batch (None, 32, 32, 128)       512       
_________________________________________________________________
conv2d_5 (Conv2D)            (None, 32, 32, 128)       147584    
_________________________________________________________________
batch_normalization_5 (Batch (None, 32, 32, 128)       512       
_________________________________________________________________
conv2d_6 (Conv2D)            (None, 32, 32, 64)        73792     
_________________________________________________________________
batch_normalization_6 (Batch (None, 32, 32, 64)        256       
_________________________________________________________________
conv2d_7 (Conv2D)            (None, 32, 32, 64)        36928     
_________________________________________________________________
batch_normalization_7 (Batch (None, 32, 32, 64)        256       
_________________________________________________________________
up_sampling2d (UpSampling2D) (None, 64, 64, 64)        0         
_________________________________________________________________
conv2d_8 (Conv2D)            (None, 64, 64, 32)        18464     
_________________________________________________________________
batch_normalization_8 (Batch (None, 64, 64, 32)        128       
_________________________________________________________________
conv2d_9 (Conv2D)            (None, 64, 64, 32)        9248      
_________________________________________________________________
batch_normalization_9 (Batch (None, 64, 64, 32)        128       
_________________________________________________________________
up_sampling2d_1 (UpSampling2 (None, 128, 128, 32)      0         
_________________________________________________________________
conv2d_10 (Conv2D)           (None, 128, 128, 1)       289       
=================================================================
Total params: 427,713
Trainable params: 426,433
Non-trainable params: 1,280

As I had mentioned earlier, I use a custom image file generator based on my 3 text files, which point to the actual image files. To build it, I combined ideas from Arindam Baidya's blog post on How to increase your small image dataset using Keras ImageDataGenerator and Nilesh's post on Writing Custom Keras Generators. Essentially the idea is to build a generator function that reads the flat files to yield batches of images, which the Keras ImageDataGenerator calls using its flow() method. To be clear, I didn't end up using any image augmentation (which is probably the biggest reason to use the ImageDataGenerator). I did try augmenting initially, but my training was not converging, which I later found to be because I was mistakenly augmenting my validation set as well. In any case, I turned off augmentation for both, since I had 10,000 images, which I figured would be more than enough to train the autoencoder. Here is the code for the generator function images_generator() and the wrapper function custom_generator().

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def images_generator(image_list):
    while True:
        indices = np.random.permutation(len(image_list))
        num_batches = len(image_list) // BATCH_SIZE
        for batch in range(num_batches):
            curr_batch = indices[batch*BATCH_SIZE:(batch + 1)*BATCH_SIZE]
            xtrain = np.empty([0, IMAGE_HEIGHT, IMAGE_WIDTH, NUM_CHANNELS],
                             dtype=np.float32)
            ytrain = np.empty([0], dtype=np.int32)
            for id_batch in curr_batch:
                image = get_source_image(IMAGE_DIR, image_list[id_batch])
                image = get_processed_image(image)
                xtrain = np.append(xtrain, [image], axis=0)
                ytrain = np.append(ytrain, [0])
            yield (xtrain, ytrain)

def custom_generator(gen, augment_images):
    if augment_images:
        keras_image_generator = tf.keras.preprocessing.image.ImageDataGenerator(
            rotation_range=15,
            width_shift_range=0.1,
            height_shift_range=0.1,
            shear_range=0.01,
            zoom_range=[0.9, 1.25],
            horizontal_flip=True,
            vertical_flip=False,
            fill_mode="reflect",
            brightness_range=[0.5, 1.5]
        )
    else:
        keras_image_generator = tf.keras.preprocessing.image.ImageDataGenerator()
    for images, labels in gen:
        kigen_batch = keras_image_generator.flow(images, labels,
                batch_size=images.shape[0])
        aug_images, labels = next(kigen_batch)
        yield (aug_images, aug_images)

In the code snippet below, the training set of 10,000 cat images is split 90/10 into a training set and validation set. We then declare custom generators on top of these image lists, and call the autoencoder's fit_generator method, training it for 50 epochs on a AWS p2.xlarge EC2 instance.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
from sklearn.model_selection import train_test_split

train_set_images = get_train_images(TRAIN_CATS)
train_images, val_images = train_test_split(train_set_images,
                                            train_size=0.9,
                                            test_size=0.1)

train_gen = images_generator(train_images)
cust_gen = custom_generator(train_gen, False)

val_gen = images_generator(val_images)
cust_val_gen = custom_generator(val_gen, False)

autoencoder = build_model_2()
print(autoencoder.summary())

checkpoint = tf.keras.callbacks.ModelCheckpoint(MODEL_FILE,
                                                save_best_only=True)
tensorboard = tf.keras.callbacks.TensorBoard(log_dir=LOG_DIR, 
                                             batch_size=BATCH_SIZE)
history = autoencoder.fit_generator(cust_gen, 
    steps_per_epoch=len(train_images) // BATCH_SIZE,
    epochs=NUM_EPOCHS,
    callbacks=[checkpoint, tensorboard],
    validation_data=cust_val_gen,
    validation_steps=len(val_images) // BATCH_SIZE)

The learning curve for the autoencoder training shown below shows that both the training and validation losses are decreasing and converging nicely.


Model Evaluation


The ModelCheckpoint callback ensures that the best model (based on the validation loss) is saved, so we can use this saved model to do our evaluations. The first evaluation I did was a few spot checks, where I rendered the input and reconstructed image (by running the input image through the trained autoencoder) side by side and computed the reconstruction error. Code for that is shown below. This code was used to produce the first two images of the cat and dog earlier in this post.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
def compute_reconstruction_error(image, image_r):
    assert(image.shape == (IMAGE_WIDTH, IMAGE_HEIGHT))
    assert(image_r.shape == (IMAGE_WIDTH, IMAGE_HEIGHT))
    return np.sqrt(np.mean(np.power(image - image_r, 2)))

test_image = os.path.join(DOG_IMAGE_DIR, "4905.jpg")
    
image_p = get_processed_image(plt.imread(test_image)

model = tf.keras.models.load_model(MODEL_FILE)

image_r = model.predict(np.expand_dims(image_p, axis=0))

image_r = np.squeeze(image_r, axis=(0, 3))
image_p = np.squeeze(image_p, axis=2)
recon_error = compute_reconstruction_error(image_p, image_r)

plt.subplot(121)
image_p = (image_p * 255).astype(np.int32)
plt.xticks([])
plt.yticks([])
plt.title("original")
plt.imshow(image_p, cmap="gray")

plt.subplot(122)
image_r = (image_r * 255).astype(np.int32)
plt.xticks([])
plt.yticks([])
plt.title("reconstructed")
plt.imshow(image_r, cmap="gray")

plt.tight_layout()
plt.show()

print("reconstruction error: {:.7f}".format(recon_error))

I followed this up by computing the reconstruction error across 2,500 cats and 2,500 dogs in the held-out test set, and plotting the results in a pair of histograms. I then computed the t-value to determine if the two distributions were different with some reasonable confidence level. The code for this is shown below. This code was used to produce the overlapping histograms earlier in the post.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
def compute_recon_errs(model, testfile, image_dir):
    num_images = 0
    images = []
    ftest = open(testfile, "r")
    for line in ftest:
        if line.startswith("#"):
            continue
        if num_images % 1000 == 0:
            print("{:d} images read".format(num_images))
        image_file = os.path.join(image_dir, line.strip())
        image = plt.imread(image_file)
        images.append(get_processed_image(image))
        num_images += 1
    
    print("{:d} images read, COMPLETE".format(num_images))
    ftest.close()
    
    print("Running model...")
    images = np.array(images)
    images_r = model.predict(images, batch_size=BATCH_SIZE)
    
    print("Computing errors...")
    errors = np.zeros((num_images,), dtype=np.float32)
    for i in range(num_images):
        image = np.squeeze(images[i], axis=2)
        image_r = np.squeeze(images_r[i], axis=2)
        errors[i] = _compute_recon_err(image, image_r)
    return errors

def compute_stats(errors):
    mu = np.mean(errors)
    sd = np.std(errors)
    return mu, sd, errors.shape[0]

def compute_t(cat_stats, dog_stats):
    cat_mu, cat_sig, cat_n = cat_stats
    dog_mu, dog_sig, dog_n = dog_stats
    mu_diff = dog_mu - cat_mu
    sd_sum = math.sqrt(((dog_sig ** 2) / dog_n) + ((cat_sig ** 2) / cat_n))
    return mu_diff / sd_sum


cat_recon_errs = compute_recon_errs(model, TEST_CATS_FILE, CAT_IMAGE_DIR)
dog_recon_errs = compute_recon_errs(model, TEST_DOGS_FILE, DOG_IMAGE_DIR)

plt.hist(cat_recon_errs, bins=50, color="b", alpha=0.3, label="cats")
plt.hist(dog_recon_errs, bins=50, color="r", alpha=0.3, label="dogs")
plt.legend(loc="best")
plt.xlabel("reconstruction error")
plt.ylabel("# of images")
plt.show()

cat_stats = compute_stats(cat_recon_errs)
dog_stats = compute_stats(dog_recon_errs)
print("cat stats:", cat_stats)
print("dog stats:", dog_stats)

t_value = compute_t(cat_stats, dog_stats)
print("t_value:", t_value)
print("lookup pvalue from: https://www.biologyforlife.com/t-test.html")

And finally, here is the snippet for computing the probability of an image being that of a cat, given the reconstruction error from the cat encoder. Essentially, I compute the value of the normal distribution at the reconstruction error point, normalized by the (maximum) value at the mean.

1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
import scipy.stats

cat_mean, cat_sd, cat_n = (0.035493888, 0.011421803, 2495)
dog_mean, dog_sd, dog_n = (0.042051364, 0.0155974245, 2494)

max_cat = scipy.stats.norm(cat_mean, cat_sd).pdf(cat_mean)

cat_recon_err = 0.04304386727035
dog_recon_err = 0.07650210412358

cat_score = scipy.stats.norm(cat_mean, cat_sd).pdf(cat_recon_error) / max_cat
dog_score = scipy.stats.norm(cat_mean, cat_sd).pdf(dog_recon_error) / max_cat

To conclude, in addition to being good at feature compression and anomaly detection, I believe autoencoders can also be useful to build complex classifiers (or feature detectors for more complex classifiers) by building on existing collections of pre-classified data. I did this experiment as a proof of concept to convince myself that this was indeed a viable approach, hopefully I have managed to convince you too.