Solved – the intuition behind the Population Stability Index

distributionsintuitionkullback-leiblerlikelihood

The "Population Stability Index" for two distributions $P$ and $Q$ is defined as the Symmetrised Kullback-Leibler divergence:

$$
\mathrm{PSI}(P,Q) = D_{KL}(P||Q) + D_{KL}(Q||P) = \sum_i(P_i-Q_i)\log\frac{P_i}{Q_i}
$$

What is the intuition behind this number?

One can always use the intuition for $D_{KL}$ and say that PSI is

the expected number of extra bits required to code samples from $P$ using a code optimized for $Q$ rather than the code optimized for $P$

plus the expected number of extra bits required to code samples from $Q$ using a code optimized for $P$ rather than the code optimized for $Q$,

but this is quite a mouthful.

Quora and UCAnalytics offer this "interpretation":

PSI < 0.1: Insignificant change (No action required)
0.1 < PSI < 0.25: Some minor change (Start worrying)
0.25 < PSI: Major shift in population (Need to delve deeper)

what is the basis for this?

Best Answer

Intuition

Kullback-Leibler Divergence can be interpreted to mean

how many bits of information we expect to lose is we use $Q$ instead of $P$.

Thus the Population Stability Index is the "roundtrip loss":

how many bits of information we expect to lose is we use $Q$ instead of $P$ and then use that again to go back to $Q$.

Values

It appears that the Population Stability Index is closely related to the G-test:

$$ \mathrm{PSI}(P,Q) = \frac{G(P,Q) + G(Q,P)}{2N} $$

(and thus can be computed using scipy.stats.power_divergence, as well as directly).

Therefore the p-values corresponding to PSI can be computed using the $\chi^2$ distribution:

import scipy.stats as st
print("            ","     ".join("DF=%d" % (df) for df in [1,2,3]))
for psi in [0.1, 0.25]:
    print "PSI=%.2f  %s" % (psi, "".join(
        " %5f" % (st.distributions.chi2.sf(psi,df)) for df in [1,2,3]))

               DF=1       DF=2       DF=3
PSI=0.10     0.751830   0.951229   0.991837
PSI=0.25     0.617075   0.882497   0.969140

Here PSI is the Population Stability Index and DF is the number of degrees of freedom ($\mathrm{DF}=n-1$ where $n$ is the number of distinct values that the variable takes).

Interestingly enough, the official "interpretation" of the PSI value completely ignores DF.

Related Solutions

Beta Distribution – Understanding the Intuition Behind Beta Distribution

The short version is that the Beta distribution can be understood as representing a distribution of probabilities, that is, it represents all the possible values of a probability when we don't know what that probability is. Here is my favorite intuitive explanation of this:

Anyone who follows baseball is familiar with batting averages—simply the number of times a player gets a base hit divided by the number of times he goes up at bat (so it's just a percentage between 0 and 1). .266 is in general considered an average batting average, while .300 is considered an excellent one.

Imagine we have a baseball player, and we want to predict what his season-long batting average will be. You might say we can just use his batting average so far- but this will be a very poor measure at the start of a season! If a player goes up to bat once and gets a single, his batting average is briefly 1.000, while if he strikes out, his batting average is 0.000. It doesn't get much better if you go up to bat five or six times- you could get a lucky streak and get an average of 1.000, or an unlucky streak and get an average of 0, neither of which are a remotely good predictor of how you will bat that season.

Why is your batting average in the first few hits not a good predictor of your eventual batting average? When a player's first at-bat is a strikeout, why does no one predict that he'll never get a hit all season? Because we're going in with prior expectations. We know that in history, most batting averages over a season have hovered between something like .215 and .360, with some extremely rare exceptions on either side. We know that if a player gets a few strikeouts in a row at the start, that might indicate he'll end up a bit worse than average, but we know he probably won't deviate from that range.

Given our batting average problem, which can be represented with a binomial distribution (a series of successes and failures), the best way to represent these prior expectations (what we in statistics just call a prior) is with the Beta distribution- it's saying, before we've seen the player take his first swing, what we roughly expect his batting average to be. The domain of the Beta distribution is (0, 1), just like a probability, so we already know we're on the right track, but the appropriateness of the Beta for this task goes far beyond that.

We expect that the player's season-long batting average will be most likely around .27, but that it could reasonably range from .21 to .35. This can be represented with a Beta distribution with parameters $\alpha=81$ and $\beta=219$:

curve(dbeta(x, 81, 219))

Beta(81, 219)

I came up with these parameters for two reasons:

The mean is $\frac{\alpha}{\alpha+\beta}=\frac{81}{81+219}=.270$
As you can see in the plot, this distribution lies almost entirely within (.2, .35)- the reasonable range for a batting average.

You asked what the x axis represents in a beta distribution density plot—here it represents his batting average. Thus notice that in this case, not only is the y-axis a probability (or more precisely a probability density), but the x-axis is as well (batting average is just a probability of a hit, after all)! The Beta distribution is representing a probability distribution of probabilities.

But here's why the Beta distribution is so appropriate. Imagine the player gets a single hit. His record for the season is now 1 hit; 1 at bat. We have to then update our probabilities- we want to shift this entire curve over just a bit to reflect our new information. While the math for proving this is a bit involved (it's shown here), the result is very simple. The new Beta distribution will be:

$\mbox{Beta}(\alpha_0+\mbox{hits}, \beta_0+\mbox{misses})$

Where $\alpha_0$ and $\beta_0$ are the parameters we started with- that is, 81 and 219. Thus, in this case, $\alpha$ has increased by 1 (his one hit), while $\beta$ has not increased at all (no misses yet). That means our new distribution is $\mbox{Beta}(81+1, 219)$, or:

curve(dbeta(x, 82, 219))

enter image description here

Notice that it has barely changed at all- the change is indeed invisible to the naked eye! (That's because one hit doesn't really mean anything).

However, the more the player hits over the course of the season, the more the curve will shift to accommodate the new evidence, and furthermore the more it will narrow based on the fact that we have more proof. Let's say halfway through the season he has been up to bat 300 times, hitting 100 out of those times. The new distribution would be $\mbox{Beta}(81+100, 219+200)$, or:

curve(dbeta(x, 81+100, 219+200))

enter image description here

Notice the curve is now both thinner and shifted to the right (higher batting average) than it used to be- we have a better sense of what the player's batting average is.

One of the most interesting outputs of this formula is the expected value of the resulting Beta distribution, which is basically your new estimate. Recall that the expected value of the Beta distribution is $\frac{\alpha}{\alpha+\beta}$. Thus, after 100 hits of 300 real at-bats, the expected value of the new Beta distribution is $\frac{81+100}{81+100+219+200}=.303$- notice that it is lower than the naive estimate of $\frac{100}{100+200}=.333$, but higher than the estimate you started the season with ($\frac{81}{81+219}=.270$). You might notice that this formula is equivalent to adding a "head start" to the number of hits and non-hits of a player- you're saying "start him off in the season with 81 hits and 219 non hits on his record").

Thus, the Beta distribution is best for representing a probabilistic distribution of probabilities: the case where we don't know what a probability is in advance, but we have some reasonable guesses.

Solved – the intuition behind conditional Gaussian distributions

Synopsis

Every statement in the question can be understood as a property of ellipses. The only property particular to the bivariate Normal distribution that is needed is the fact that in a standard bivariate Normal distribution of $X,Y$--for which $X$ and $Y$ are uncorrelated--the conditional variance of $Y$ does not depend on $X$. (This in turn is an immediate consequence of the fact that lack of correlation implies independence for jointly Normal variables.)

The following analysis shows precisely what property of ellipses is involved and derives all the equations of the question using elementary ideas and the simplest possible arithmetic, in a way intended to be easily remembered.

Circularly symmetric distributions

The distribution of the question is a member of the family of bivariate Normal distributions. They are all derived from a basic member, the standard bivariate Normal, which describes two uncorrelated standard Normal distributions (forming its two coordinates).

Figure 1: the standard bivariate normal distribution

The left side is a relief plot of the standard bivariate normal density. The right side shows the same in pseudo-3D, with the front part sliced away.

This is an example of a circularly symmetric distribution: the density varies with distance from a central point but not with the direction away from that point. Thus, the contours of its graph (at the right) are circles.

Most other bivariate Normal distributions are not circularly symmetric, however: their cross-sections are ellipses. These ellipses model the characteristic shape of many bivariate point clouds.

Figure 2: another bivariate normal distribution, plotted

These are portraits of the bivariate Normal distribution with covariance matrix $\Sigma = \left(\begin{array}{cc} 1 & -\frac{2}{3} \\ -\frac{2}{3} & 1 \\\end{array}\right).$ It is a model for data with correlation coefficient $-2/3$.

How to Create Ellipses

An ellipse--according to its oldest definition--is a conic section, which is a circle distorted by a projection onto another plane. By considering the nature of projection, just as visual artists do, we may decompose it into a sequence of distortions that are easy to understand and calculate with.

First, stretch (or, if necessary, squeeze) the circle along what will become the long axis of the ellipse until it is the correct length:

Step 1: stretch

Next, squeeze (or stretch) this ellipse along its minor axis:

Step 2: squeeze

Third, rotate it around its center into its final orientation:

Step 3: rotate

Finally, shift it to the desired location:

Step 4: shift

These are all affine transformations. (In fact, the first three are linear transformations; the final shift makes it affine.) Because a composition of affine transformations is (by definition) still affine, the net distortion from the circle to the final ellipse is an affine transformation. But it can be somewhat complicated:

Composite transformation

Notice what happened to the ellipse's (natural) axes: after they were created by the shift and squeeze, they (of course) rotated and shifted along with the axis itself. We easily see these axes even when they are not drawn, because they are axes of symmetry of the ellipse itself.

We would like to apply our understanding of ellipses to understanding distorted circularly symmetric distributions, like the bivariate Normal family. Unfortunately, there is a problem with these distortions: they do not respect the distinction between the $x$ and $y$ axes. The rotation at step 3 ruins that. Look at the faint coordinate grids in the backgrounds: these show what happens to a grid (of mesh $1/2$ in both directions) when it is distorted. In the first image the spacing between the original vertical lines (shown solid) is doubled. In the second image the spacing between the original horizontal lines (shown dashed) is shrunk by a third. In the third image the grid spacings are not changed, but all the lines are rotated. They shift up and to the right in the fourth image. The final image, showing the net result, displays this stretched, squeezed, rotated, shifted grid. The original solid lines of constant $x$ coordinate no longer are vertical.

The key idea--one might venture to say it is the crux of regression--is that there is a way in which the circle can be distorted into an ellipse without rotating the vertical lines. Because the rotation was the culprit, let's cut to the chase and show how to created a rotated ellipse without actually appearing to rotate anything!

Skewed ellipse

This is a skew transformation. It does two things at once:

It squeezes in the $y$ direction (by an amount $\lambda$, say). This leaves the $x$-axis alone.
It lifts any resulting point $(x,\lambda y)$ by an amount directly proportional to $x$. Writing that constant of proportionality as $\rho$, this sends $(x,\lambda y)$ to $(x, \lambda y+\rho x)$.

The second step lifts the $x$-axis into the line $y=\rho x$, shown in the previous figure. As shown in that figure, I want to work with a special skew transformation, one that effectively rotates the ellipse by 45 degrees and inscribes it into the unit square. The major axis of this ellipse is the line $y=x$. It is visually evident that $|\rho| \le 1$. (Negative values of $\rho$ tilt the ellipse down to the right rather than up to the right.) This is the geometric explanation of "regression to the mean."

Choosing an angle of 45 degrees makes the ellipse symmetric around the square's diagonal (part of the line $y=x$). To figure out the parameters of this skew transformation, observe:

The lifting by $\rho x$ moves the point $(1,0)$ to $(1,\rho)$.
The symmetry around the main diagonal then implies the point $(\rho, 1)$ also lies on the ellipse.

Where did this point start out?

The original (upper) point on the unit circle (having implicit equation $x^2+y^2=1$) with $x$ coordinate $\rho$ was $(\rho, \sqrt{1-\rho^2})$.
Any point of the form $(\rho, y)$ first got squeezed to $(\rho, \lambda y)$ and then lifted to $(\rho, \lambda y + \rho\times\rho)$.

The unique solution to the equation $(\rho, \lambda \sqrt{1-\rho^2} + \rho^2) = (\rho, 1)$ is $\lambda = \sqrt{1-\rho^2}$. That is the amount by which all distances in the vertical direction must be squeezed in order to create an ellipse at a 45 degree angle when it is skewed vertically by $\rho$.

To firm up these ideas, here is a tableau showing how a circularly symmetric distribution is distorted into distributions with elliptical contours by means of these skew transformations. The panels show values of $\rho$ equal to $0,$ $3/10,$ $6/10,$ and $9/10,$ from left to right.

Tableau

The leftmost figure shows a set of starting points around one of the circular contours as well as part of the horizontal axis. Subsequent figures use arrows to show how those points are moved. The image of the horizontal axis appears as a slanted line segment (with slope $\rho$). (The colors represent different amounts of density in the different figures.)

Application

We are ready to do regression. A standard, elegant (yet simple) method to perform regression is first to express the original variables in new units of measurement: we center them at their means and use their standard deviations as the units. This moves the center of the distribution to the origin and makes all its elliptical contours slant 45 degrees (up or down).

When these standardized data form a circular point cloud, the regression is easy: the means conditional on $x$ are all $0$, forming a line passing through the origin. (Circular symmetry implies symmetry with respect to the $x$ axis, showing that all conditional distributions are symmetric, whence they have $0$ means.) As we have seen, we may view the standardized distribution as arising from this basic simple situation in two steps: first, all the (standardized) $y$ values are multiplied by $\sqrt{1-\rho^2}$ for some value of $\rho$; next, all values with $x$-coordinates are vertically skewed by $\rho x$. What did these distortions do to the regression line (which plots the conditional means against $x$)?

The shrinking of $y$ coordinates multiplied all vertical deviations by a constant. This merely changed the vertical scale and left all conditional means unaltered at $0$.
The vertical skew transformation added $\rho x$ to all conditional values at $x$, thereby adding $\rho x$ to their conditional mean: the curve $y=\rho x$ is the regression curve, which turns out to be a line.

Similarly, we may verify that because the $x$-axis is the least squares fit to the circularly symmetric distribution, the least squares fit to the transformed distribution also is the line $y=\rho x$: the least-squares line coincides with the regression line.

These beautiful results are a consequence of the fact that the vertical skew transformation does not change any $x$ coordinates.

We can easily say more:

The first bullet (about shrinking) shows that when $(X,Y)$ has any circularly symmetric distribution, the conditional variance of $Y|X$ was multiplied by $\left(\sqrt{1-\rho^2}\right)^2 = 1 - \rho^2$.
More generally: the vertical skew transformation rescales each conditional distribution by $\sqrt{1-\rho^2}$ and then recenters it by $\rho x$.

For the standard bivariate Normal distribution, the conditional variance is a constant (equal to $1$), independent of $x$. We immediately conclude that after applying this skew transformation, the conditional variance of the vertical deviations is still a constant and equals $1-\rho^2$. Because the conditional distributions of a bivariate Normal are themselves Normal, now that we know their means and variances, we have full information about them.

Finally, we need to relate $\rho$ to the original covariance matrix $\Sigma$. For this, recall that the (nicest) definition of the correlation coefficient between two standardized variables $X$ and $Y$ is the expectation of their product $XY$. (The correlation of $X$ and $Y$ is simply declared to be the correlation of their standardized versions.) Therefore, when $(X,Y)$ follows any circularly symmetric distribution and we apply the skew transformation to the variables, we may write

$$\varepsilon = Y - \rho X$$

for the vertical deviations from the regression line and notice that $\varepsilon$ must have a symmetric distribution around $0$. Why? Because before the skew transformation was applied, $Y$ had a symmetric distribution around $0$ and then we (a) squeezed it and (b) lifted it by $\rho X$. The former did not change its symmetry while the latter recentered it at $\rho X$, QED. The next figure illustrates this.

3D plot showing conditional distributions and the least-squares line

The black lines trace out heights proportional to the conditional densities at various regularly-spaced values of $x$. The thick white line is the regression line, which passes through the center of symmetry of each conditional curve. This plot shows the case $\rho = -1/2$ in standardized coordinates.

Consequently

$$\mathbb{E}(XY) = \mathbb{E}\left(X(\rho X + \varepsilon)\right) = \rho\mathbb{E}(X^2) + \mathbb{E}(X\varepsilon) = \rho(1) + 0=\rho.$$

The final equality is due to two facts: (1) because $X$ has been standardized, the expectation of its square is its standardized variance, equal to $1$ by construction; and (2) the expectation of $X\varepsilon$ equals the expectation of $X(-\varepsilon)$ by virtue of the symmetry of $\varepsilon$. Because the latter is the negative of the former, both must equal $0$: this term drops out.

We have identified the parameter of the skew transformation, $\rho$, as being the correlation coefficient of $X$ and $Y$.

Conclusions

By observing that any ellipse may be produced by distorting a circle with a vertical skew transformation that preserves the $x$ coordinate, we have arrived at an understanding of the contours of any distribution of random variables $(X,Y)$ that is obtained from a circularly symmetric one by means of stretches, squeezes, rotations, and shifts (that is, any affine transformation). By re-expressing the results in terms of the original units of $x$ and $y$--which amount to adding back their means, $\mu_x$ and $\mu_y$, after multiplying by their standard deviations $\sigma_x$ and $\sigma_y$--we find that:

The least-squares line and the regression curve both pass through the origin of the standardized variables, which corresponds to the "point of averages" $(\mu_x,\mu_y)$ in original coordinates.
The regression curve, which is defined to be the locus of conditional means, $\{(x, \rho x)\},$ coincides with the least-squares line.
The slope of the regression line in standardized coordinates is the correlation coefficient $\rho$; in the original units it therefore equals $\sigma_y \rho / \sigma_x$.

Consequently the equation of the regression line is

$$y = \frac{\sigma_y\rho}{\sigma_x}\left(x - \mu_x\right) + \mu_y.$$

The conditional variance of $Y|X$ is $\sigma_y^2(1-\rho^2)$ times the conditional variance of $Y'|X'$ where $(X',Y')$ has a standard distribution (circularly symmetric with unit variances in both coordinates), $X'=(X-\mu_X)/\sigma_x$, and $Y'=(Y-\mu_Y)/\sigma_Y$.

None of these results is a particular property of bivariate Normal distributions! For the bivariate Normal family, the conditional variance of $Y'|X'$ is constant (and equal to $1$): this fact makes that family particularly simple to work with. In particular:

Because in the covariance matrix $\Sigma$ the coefficients are $\sigma_{11}=\sigma_x^2,$ $\sigma_{12}=\sigma_{21}=\rho\sigma_x\sigma_y,$ and $\sigma_{22}=\sigma_y^2,$ the conditional variance of $Y|X$ for a bivariate Normal distribution is

$$\sigma_y^2(1-\rho^2)=\sigma_{22}\left(1-\left(\frac{\sigma_{12}}{\sqrt{\sigma_{11}\sigma_{22}}}\right)^2\right)=\sigma_{22} - \frac{\sigma_{12}^2}{\sigma_{11}}.$$

Technical Notes

The key idea can be stated in terms of matrices describing the linear transformations. It comes down to finding a suitable "square root" of the correlation matrix for which $y$ is an eigenvector. Thus:

$$\left(\begin{array}{cc} 1 & \rho \\ \rho & 1 \\\end{array}\right) = \mathbb{A}\mathbb{A}'$$

where

$$\mathbb{A} = \left(\begin{array}{cc} 1 & 0 \\ \rho & \sqrt{1-\rho^2} \\\end{array}\right).$$

A much better known square root is the one initially described (involving a rotation instead of a skew transformation); it is the one produced by a singular value decomposition and it plays a prominent role in principal components analysis (PCA):

$$\left(\begin{array}{cc} 1 & \rho \\ \rho & 1 \\\end{array}\right) = \mathbb{B}\mathbb{B}';$$

$$\mathbb{B} = \mathbb{Q} \left( \begin{array}{cc} \sqrt{\rho +1} & 0 \\ 0 & \sqrt{1-\rho } \\ \end{array} \right)\mathbb{Q}'$$

where $\mathbb{Q} = \left( \begin{array}{cc} \frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\ \end{array} \right)$ is the rotation matrix for a $45$ degree rotation.

Thus, the distinction between PCA and regression comes down to the difference between two special square roots of the correlation matrix.