Histograms and kernel density estimation with Rapaio toolbox

Histogram and Kernel Density Estimation Tutorial

Data set used

First we need to load a frame with data. For convenience we don't use input/output facilities. Instead, we load a set of data already built-in the library.

We will use the classical Pearson father-son data set.

>>summary("pearson", [Father, Son])
rows: 1078, cols: 2
          Father               Son 
   Min. : 59.000     Min. : 58.500 
1st Qu. : 65.800  1st Qu. : 66.900 
 Median : 67.800   Median : 68.600 
   Mean : 67.687     Mean : 68.684 
2nd Qu. : 69.600  2nd Qu. : 70.500 
   Max. : 75.400     Max. : 78.400 
                                   

We have two continuous random variables which observes height of pairs of father and sons.

The 5-number summaries gives us some hints on how the distribution of the values for the 2 random variables looks like. We note that the mean and median are very close which leads us to think that the distribution of the values are somehow symmetric.

Histograms

One of the most important and usual tools to obtain information about the distribution of a random variable is histogram.

In statistics, a histogram is a graphical representation of the distribution of data. It is an estimate of the probability distribution of a continuous variable and was first introduced by Karl Pearson. A histogram is a representation of tabulated frequencies, shown as adjacent rectangles, erected over discrete intervals (bins), with an area equal to the frequency of the observations in the interval.

We draw histograms with Rapaio toolbox in the following way:

        draw(new Histogram(df.getCol("Father")));

The height of a rectangle is also equal to the frequency density of the interval, i.e., the frequency divided by the width of the interval. The total area of the histogram is equal to the number of data.

The previous code uses the simplest form of building a histogram. There are some parameters which can help somebody to tune the output.

One can change the number of bins.

        draw(new Histogram(df.getCol("Father"), 100, false));

Note that on the vertical axis we found the count of the elements which are held by the bins that are displayed. We can change how the heights of the bins are computed into densities which makes the total area under curve to be 1. That feature is a key property of a probability density function, also.

draw(new Histogram(df.getCol("Father"), 30, true));

The histogram is useful but have a weak point. Its weak point lies into it's flexibility given by the number of bins. There is no "best" number of bins, and different bin sizes can reveal different features of the data.

Some theoreticians have attempted to determine an optimal number of bins, but these methods generally make strong assumptions about the shape of the distribution. Depending on the actual data distribution and the goals of the analysis, different bin widths may be appropriate, so experimentation is usually needed to determine an appropriate width.

An alternative to the histogram is kernel density estimation, which uses a kernel to smooth samples. This will construct a smooth probability density function, which will in general more accurately reflect the underlying variable.

One can draw also the kernel density approximation, over a histogram or as a separate plot.

        Vector col = df.getCol("Father");
        Plot plot = new Plot();
        HistogramBars hb = new HistogramBars(plot, col);
        hb.opt().setColorIndex(new IndexVector("rainbow", 1, 255, 1));
        plot.add(hb);
        plot.add(new DensityLine(plot, col));
        draw(plot);

In statistics, kernel density estimation (KDE) is a non-parametric way to estimate the probability density function of a random variable. Kernel density estimation is a fundamental data smoothing problem where inferences about the population are made, based on a finite data sample. In some fields such as signal processing and econometrics it is also termed the Parzen–Rosenblatt window method, after Emanuel Parzen and Murray Rosenblatt, who are usually credited with independently creating it in its current form.

In it's default implementation, used without parameters, the Rapaio toolbox build a kernel density estimation with Gaussian kernels and with bandwidth approximated by Silverman's rule of thumb.

However one can use a different value for bandwidth in order to obtain a smooth or less smooth approximation of the density function.

Another thing one can try with kernel density estimator is to change the kernel function, which is the function used to disperse in a small neighborhood around a point the probability mass initially assigned to the points of the discrete sample. The kernel choices are: uniform, triangular, Epanechnikov, biweight, triweight, tricube, Gaussian and cosine. Of course, it is easy to implement your own smoothing method once you implement a custom kernel function.

We could agree that my implementation of kernel function is ugly and maybe no so useful, however you have to know that the purpose of Rapaio toolbox is to give you sharp and precise standard tools and, in the same time, the opportunity to experiment with your owns.

One final graph will show the kernel approximation of density functions for both numerical variables that we have: father's heights and son's heights.

Blue line represents density approximation of father's heights, red line represents density approximation of son's heights.

Note: the sole purpose of this tutorial is to show what and how it can be done with Rapaio toolbox library.

>>>This tutorial is generated with Rapaio document printer facilities.<<<

Pearson product-moment and Spearman's rank correlations with Rapaio

Introduction

This tutorial presents you the correlations tools offered by Rapaio library.
We will use the classical iris data set. The numerical columns of this dataset are:

Frame df = Datasets.loadIrisDataset();
df = ColFilters.retainNumeric(df);
names(df);

>>names("iris")
sepal-length
sepal-width
petal-length
petal-width

Pearson product-moment correlation

Pearson product-moment correlation measures the linear correlation between two random variables. Among other type of correlation measures, the Pearson product-moment detects only linear correlations.

Definition

Pearson product-moment coefficient measures the linear correlation between two random variables \(X\) and \(Y\), giving a value between +1 and −1 inclusive, where 1 is total positive correlation, 0 is no correlation, and −1 is negative correlation.
Pearson's correlation coefficient when applied to a population is commonly represented by the Greek letter \( ho\) (rho) and may be referred to as the population correlation coefficient or the population Pearson correlation coefficient. The formula for \(\rho\) is:
$$ \rho_{X,Y} = \frac{cov(X,Y)}{\sigma_X\sigma_Y} = \frac{E[(X-\mu_X)(Y-\mu_Y)]}{\sigma_X\sigma_Y} $$Pearson's correlation coefficient when applied to a sample is commonly represented by the letter \(r\) and may be referred to as the sample correlation coefficient or the sample Pearson correlation coefficient. We can obtain a formula for \(r\) by substituting estimates of the covariances and variances based on a sample into the formula above. That formula for \(r\) is:
$$ r = \frac{\sum ^n _{i=1}(X_i-\bar{X})(Y_i-\bar{Y})}{\sqrt{\sum ^n _{i=1}(X_i - \bar{X})^2} \sqrt{\sum ^n _{i=1}(Y_i - \bar{Y})^2}} $$The interpretation of a correlation coefficient depends on the context and purposes. A correlation of 0.8 may be very low if one is verifying a physical law using high-quality instruments, but may be regarded as very high in the social sciences where there may be a greater contribution from complicating factors.

Usa Rapaio for Pearson correlation

Rapaio library allows one to compute Pearson \(r\) for more then one vector at a time. Thus the result will be a matrix with computed \(r\) values between vectors,using vectors index position as indexes in resulted matrix.

        PearsonRCorrelation corr = new PearsonRCorrelation(df);
        summary(corr);

pearson[[sepal-length, sepal-width, petal-length, petal-width]] - Pearson product-moment correlation coefficient
   1.sepal-length 2.sepal-width 3.petal-length 4.petal-width 
1.              x     -0.109369       0.871754      0.817954 
2.      -0.109369             x      -0.420516     -0.356544 
3.       0.871754     -0.420516              x      0.962757 
4.       0.817954     -0.356544       0.962757             x 

We can spot with eas that many of the attributes are linearly correlated. As a sample we find from the correlation summary that petal-length and petal-width have a very strong linear correlation. Let's check this intuition with a plot:

Another \(r\) coefficient which have a value close to \(1\) is between sepal-length and petal-length. Let's check that with a plot, also:

Finally, we plot again, but this time using a coefficient which is closer to 0, which could mean that the variables are not linearly correlated. Such a value for correlation we have between sepal-length and sepal-width.

Spearman's rank correlation coefficient

often denoted by the Greek letter \(\rho\) (rho) or as \(r_s\), is a nonparametric measure of statistical dependence between two variables. It assesses how well the relationship between two variables can be described using a monotonic function. If there are no repeated data values, a perfect Spearman correlation of \(+1\) or \(−1\) occurs when each of the variables is a perfect monotone function of the other.

Definition

The Spearman correlation coefficient is defined as the Pearson correlation coefficient between the ranked variables. For a sample of size \(n\), the \(n\) raw scores \(X_i\), \(Y_i\) are converted to ranks \(x_i\), \(y_i\), and \(\rho\) is computed from these:
$$ \rho = \frac{\sum_i(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum_i (x_i-\bar{x})^2 \sum_i(y_i-\bar{y})^2}} $$Identical values (rank ties or value duplicates) are assigned a rank equal to the average of their positions in the ascending order of the values.

Use Rapaio to compute Spearman's rank correlation

Rapaio library allows one to compute Spearman \(\rho\) for more then one vector at a time. Thus the result will be a matrix with computed \(\rho\) values between vectors,using vectors index position as indexes in resulted matrix.

spearman[[sepal-length, sepal-width, petal-length, petal-width]] - Spearman's rank correlation coefficient
   1.sepal-length 2.sepal-width 3.petal-length 4.petal-width 
1.              x     -0.159457       0.881386      0.834421 
2.      -0.159457             x      -0.303421     -0.277511 
3.       0.881386     -0.303421              x      0.936003 
4.       0.834421     -0.277511       0.936003             x 

pearson[[sepal-length, sepal-width, petal-length, petal-width]] - Pearson product-moment correlation coefficient
   1.sepal-length 2.sepal-width 3.petal-length 4.petal-width 
1.              x     -0.109369       0.871754      0.817954 
2.      -0.109369             x      -0.420516     -0.356544 
3.       0.871754     -0.420516              x      0.962757 
4.       0.817954     -0.356544       0.962757             x 

Rapaio - statistical, data mining and machine learning toolbox

A year ago I participated to a machine learning contest. The contest was about predicting what movies would like each user. The results of the contest are not important. But due to the fact that the contest was held for one year period, I had the chance to try a lot of ideas along the way. I knew nothing about machine learning, and statistics, and I tried somehow to see what an engineer ca do without prior specific knowledge.

What I have earned was that I discovered a whole new world where many precise tools are available, and almost all of them can be useful only if their results are properly interpreted. I discovered a world were, at least that's what I think now, there is no chance to do something without a proper human judgment. And what I found astonishing is that all those available tools are no more than precise common sense. I loved it. I feel like I am captured in this new world and in the same time, completely free to explore, find and understand how the real world works. I started to buy books and study them line by line.

As usual, when I find something worth to learn, I always choose the hard way. And for me, the hard way in this case means to try to implements all those tools in order to understand better their limitations, theirs reasoning and power. The result of this starting efforts are now materialized into a project, an open source project.

Rapaio library is a statistics, data mining and machine learning toolbox. Now it is work-in-progress and most probably will stay this way for a long period of time. That does not mean it will not be usable. It means that ahead is a myriad of tools and knowledge waiting to be implemented and used.

The source code is available on githib.
https://github.com/padreati/rapaio


And to have a preview on what this library is doing I started to create some tutorials on various topics, where I present how this library can help you. Those tutorials are strongly connected with the source code itself. The code which generates those tutorial are part of the library itself and that will have among other benefits, the fact that those tutorials will be up-to-date and ready to use in practice.

Rapaio Tutorials are availableon github also.
http://padreati.github.io/rapaio/

Why variance estimator is biased when uses n as denominator.

Given  \( X_1, .. , X_n \) a set of independent and identically distributed random variables with fixed, existent and unknown mean \( \mu \) and variance \(\sigma^2\), one usual problem which needs to be solved is to estimate population variance \(\sigma^2\). One well known estimators for that, which we will denote as \(S_n^2\) is defined by the following formula.
$$ S_n^2 = E\left[\frac{1}{n}\sum_{i=1}^{n}(X_i-\bar{X_n})^2\right] $$
where \(\bar{X_n} = \frac{1}{n}\sum_{i=1}^{n}X_i \) is the population mean estimator.

It is a well know fact that the estimator for population variance \(S_n^2\) is biased when it uses \(n\) as denominator. By biased, we mean that it's expected value does not equals with the theoretical value of population variance \(\sigma^2\). This fact is stated in almost all statistics books and this is the reason why this estimator is not considered THE population variance estimator. Even so, I saw no complete/clear proof of this fact. On wikipedia site there are two proofs but none convinced me enough to accept that and this problem bothered me for quite a few weeks.

Finally I arrived at a simpler (in my humble opinion) proof of that fact: this population variance estimator is biased.

To proceed with the proof we have to assert that for a given random variable \(X\), the expected value of \(X^2\) equals squared mean plus variation.
$$ EX^2 = E(X-\mu+\mu)^2 = E(X-\mu)^2 + 2\mu E(X-\mu) + E\mu^2 = \sigma^2+\mu^2 $$
Another needed assertion is the fact that expected value of the product of two independent random variables equals the product of the expected values of those variables. $$ EXY = EXEY $$.

Now we proceed to develop the definition of \(S_n^2\)

$$ S_n^2 = E\left[\frac{1}{n}\sum_{i=1}^{n}(X_i - \bar{X_n})^2\right] = E\left[\frac{1}{n}\sum_{i=1}^{n}(X_i^2 - 2X_i\bar{X_n} + \bar{X_n}^2)\right] = \frac{1}{n}\sum_{i=1}^{n}EX_i^2 - \frac{2}{n}\sum_{i=1}^{n}EX_i\bar{X_n} + \frac{1}{n}\sum_{i=1}^{n}E\bar{X_n}^2 $$

Using preliminary results we have

$$ = \mu^2 + \sigma^2 - \frac{2}{n}\sum_{i=1}^{n}EX_i\left(\frac{1}{n}\sum_{j=1}^{n}X_j\right) + \mu^2 + \frac{\sigma^2}{n} = 2\mu^2 + \sigma^2 + \frac{\sigma^2}{n} - \frac{2}{n^2}\sum_{i=1}^{n}\sum_{j=1}^{n}EX_iX_j $$

In the last equation we study the sum of sums of products \(\sum_{i=1}^{n}\sum_{j=1}^{n}EX_iX_j\). This expression has \(n^2\) terms and, among these terms, there are \(n\) terms with equal indexes for which we empower first preliminary result. Also, the same expression contains \(n^2-n\) terms with different indexes, where we can empower second preliminary. Thus we can develop further the equation of \(S_n^2\)

$$ S_n^2 = 2\mu^2 + \sigma^2 + \frac{\sigma^2}{n} - \frac{2}{n^2}\left[n(\mu^2+\sigma^2) + (n^2-n)\mu^2\right]
= 2\mu^2 + \sigma^2 + \frac{\sigma^2}{n} - \frac{2}{n^2}\left[n\mu^2+n\sigma^2 + n^2\mu^2 - n\mu^2\right] $$
$$ S_n^2 = 2\mu^2 + \sigma^2 + \frac{\sigma^2}{n} - 2\mu^2 - \frac{2\sigma^2}{n} = \sigma^2 - \frac{\sigma^2}{n} $$

Finally, we find that the estimator is biased, i.e. its expected value is not equal with theoretical population variance \(\sigma^2\).
$$ S_n^2 = \sigma^2\frac{n-1}{n} \leq \sigma^2 $$