There are many ways to measure the similarity between observations. The most appropriate measure has to be intuitive and avoid computational complexity. This section reviews three similarity measures:

The Manhattan distance is defined by the absolute distance between two variables or vectors, {x_i} and {y_i}, of the same size:

The implementation is generic enough to compute the distance between two arrays of elements of different types as long as an implicit conversion between each of these types to Double values is already defined, as shown here:

def manhattan[T <% Double, U <% Double](x: Array[T], y: Array[U]): Double = (x, y).zipped.foldLeft(0.0)((s, t) => s + Math.abs(t._1 - t._2))

The ubiquitous Euclidean distance is defined as the square of the distance between two vectors, {x_i} and {y_i}, of the same size:

def euclidean[T <% Double, U <% Double](x: Array[T], y: Array[U]): Double =  Math.sqrt((x, y).zipped.foldLeft(0.0)((s, t) => { val d = t._1 - t._2; s + d*d} ))

The cosine distance is defined as the cosine of an angle between two vectors, {x_i} and {y_i}, of the same size:

In this implementation, the computation of the dot product and the norms for each dataset is done simultaneously using the tuple within the fold method:

def cosine[T <% Double, U <% Double](x: Array[T], y: Array[U]): Double = {
  val zeros = (0.0, 0.0, 0.0)
  val norms = (x, y).zipped.foldLeft(zeros)((s, t) => 
     (s._1 + t._1*t._2, s._2 + t._1*t._1, s._3 + t._2*t._2))
  norms._1/Math.sqrt(norms._2*norms._3)
}

def dot (x:Array[Double], y:Array[Double]): Array[Double] = 
    x.zip(y).map( x => f(x._1, x._2) )
An functional alternative is to use the Tuple2.zipped method.
def dot(x:Array[Double], y:Array[Double]): Array[Double] = (x, y).zipped map ( _ * _)

The main advantage of the K-means algorithm (and the reason for its popularity) is its simplicity [4:3].

Note

Let's consider K clusters {C_k} with means {m_k}. The K-means algorithm is indeed an optimization problem, the objective of which is to minimize the reconstruction or the total error defined as the total sum of distance.

The steps of the iterative algorithm are:

The configuration of the K clusters consists of defining the following parameters for the K-means algorithm: number of K clusters, the distance metrics, the maximum number of iterations, and the initial value of the cluster's centroid.

class Cluster[T <% Double](val center: DblVector) {
  val members = new ListBuffer[Int]

object Cluster {
  def apply[T <% Double](c:DblVector):Cluster[T] = new Cluster[T](c)
}

def += (n:Int): Unit = members.append(n)
def moveCenter(xt: XTSeries[Array[T]): Cluster[T] ={
   val sums = members.map(xt(_).map(_.toDouble)).toList
                     .transpose
                     .map( _.sum)
  Cluster[T](sums.map( _ / members.size).toArray)
}

def stdDev(xt: XTSeries[Array[T]], distance: (DblVector, Array[T]) => Double): Double = {
   Stats[Double](members.map(xt( _))
             .map( distance(center, _)).toArray).stdDev
}

class KMeans[T <% Double](K: Int, maxIters: Int, distance: (DblVector,Array[T]) => Double)(implicit order: Ordering[T], m: Manifest[T]) extends PipeOperator[XTSeries[Array[T]], List[Cluster[T]]] {
  def |> : PartialFunction[XTSeries[Array[T]], List[Cluster[T]]]
  def initialize(xt:XTSeries[Array[T]]): List[Cluster[T]]

  def initialize(xt:XTSeries[Array[T]]): List[Cluster[T]]={
   val stats = statistics(xt) //1
   val maxSDevDim = Range(0,stats.size).maxBy (stats( _ ).stdDev)//2
   val rankedObs = xt.zipWithIndex
                     .map(x=> (x._1(maxSDevDim), x._2)) //2
                     .sortWith( _._1  < _._1) //3
   val halfSegSize = ((rankedObs.size>>1)/K).floor.toInt //4
   val centroids = rankedObs.filter(isContained( _, halfSegSize, rankedObs.size) ).map(n => xt(n._2)) //6
   Range(0, K).foldLeft(List[Cluster[T]]())((xs, i) => Cluster[T](centroids(i)) :: xs) //7
}

def isContained(t: (T,Int), hSz: Int, dim: Int): Boolean = 
    ((t._2 % hSz == 0) && (t._2 %(hSz<<1) != 0)

The second step in the K-means algorithm is the assignment of the observations to the clusters for which the centroids have been initialized in step 1. This feat is accomplished by the private assignToClusters method:

def assignToClusters(xt: XTSeries[Array[T]], clusters: List[Cluster[T]], membership: Array[Int]): Int = {
  xt.toArray
    .zipWithIndex
    .filter(x => { //1
       val nearestCluster = getNearestCluster(clusters, x._1)//2
       val reassigned = nearestCluster != membership(x._2) 
       clusters(nearestCluster) += x._2 //3
       membership(x._2) = nearestCluster //4
       reassigned
     }).size
}

The core of the assignment of observations to each cluster is the filter on the time series (line 1). The filter computes the index of the closest cluster and checks whether the observation is to be reassigned (line 2). The observation at the index x._2 is added to the nearest cluster, clusters(nearestCluster) (line 3). The current membership of the observations is then updated (line 4).

The cluster closest to an observation data is computed by the getNearestCluster method as follows:

def getNearestCluster(clusters: List[Cluster[T]], x:Array[T]): Int={
  clusters.zipWithIndex..foldLeft((Double.MaxValue,0))((p,c) => { 
      val measure = distance(c._1.center, x)
      if(measure < p._1) (measure, c._2) else p
   })._2

A fold is used to extract from the list of clusters the cluster that is closest to the observation x using the distance metric defined in the K-means constructor.

The final step is to implement the iterative computation of the reconstruction error. In this implementation, the iteration terminates when no more observations are reassigned to different clusters. As with other data processing units, the extraction of K-means clusters is encapsulated by the pipe operator |>, so that clustering can be integrated into a workflow using dependency injection described in the Dependency Injection section in Chapter 2, Hello World!.

The generation of the K clusters is executed by the data transformation |>:

def |> :PartialFunction[XTSeries[Array[T]], List[Cluster[T]]] = {
  case xt: XTSeries[Array[T]] if(xt.size>2 && xt(0).size>0) => {
    val clusters = initialize(xt)  //1

    if( clusters.isEmpty) List.empty
    else  {
      val membership = Array.fill(xt.size)(0)
      val reassigned = assignToClusters(xt,clusters,membership)//2
      var newClusters: List[Cluster[T]] = List.empty
      Range(0, maxIters).find( _ => {
        newClusters = clusters.map( c => {
          if( c.size > 0) c.moveCenter(xt, dimension(xt)) 
          else clusters.filter( _.size > 0)
                       .maxBy( _.stdDev(xt, distance))
        }) //3
        assignToClusters(xt, newClusters, membership) == 0
      }) match {
        case Some(index) => newClusters
        case None => { … }
    } //4
  }
}

As described in the algorithm overview section, the main method initializes the membership for all the observations (line 1), creates and initializes the clusters, and assigns the observations to clusters using the assignToClusters method (line 2). The iteration updates the content of each cluster using the moveCenter method, by assigning new observations to the cluster with the highest standard deviation (line 3). The iterative loop exits when no more reassignment is needed (line 4).

The companion object for KMeans implements the apply constructor and the computation of the stdDev standard deviation for each cluster. The default constructor uses the Euclidean distance:

def apply[T <% Double](K: Int, maxIters: Int)(implicit order: Ordering[T], m: Manifest[T]): KMeans[T] = new KMeans[T](K, maxIters, euclidean)
def stdDev[T](c: List[Cluster[T]], xt: XTSeries[Array[T]]): List[Double] =  c.map( _.stdDev(xt))

The stdDev method computes the standard deviation of the distances between each data point that belongs to a c cluster and its centroid.

Note that ordering a trait and Manifest have to be provided in the apply constructor because there is no guarantee that such capabilities are provided in runtime by the client code.

A model with a significant number of features (high dimensions) requires a larger number of observations in order to extract robust clusters. K-means clustering with very small datasets, of size less than 50, produces models with high bias and a limited number of clusters that are affected by the order of observations [4:5]. I have been using the following simple empirical rule of thumb for a training set of size n, expected K clusters, and N features: n < K.N.

Whichever empirical rule you follow, such a restriction is particularly an issue for analyzing stocks using historical quotes. Let's consider our examples of using technical analysis to categorize stocks according to their price behavior over a period of 1 year (or approximately 250 trading days). The dimension of the problem is 250 (250 daily closing prices). The number of stocks (observations) would have exceeded several hundred!

Price model for K-means clustering

There are options to get around this limitation and shrink the number of observations; among them are:

These approaches are workaround solutions at best, used for the sake of this tutorial, but they are not recommended for actual commercial analytical applications. The principal component analysis introduced in the last section of this chapter is one of the most reliable dimension reduction techniques.

The objective is to extract clusters from a set of stock price actions during a period of time between January 1 and Dec 31, 2013 as features. For this test, 127 stocks are randomly selected from the S&P 500 list. The following chart visualizes the behavior of the normalized price of a subset of these 127 stocks:

Price action of stocks used in K-means clustering

The key is to select the appropriate features prior to clustering and the time window to operate on. It would make sense to consider the entire historical price over the 252 trading days as a feature. However, the number of observations (stocks) is too limited to use the entire price range. The (SAMPLES = 50) observations are the stock closing price for each trading session between the 80^th and 130^th days. The adjusted daily closing prices are normalized using the minimum and maximum values.

First, let's create a simple function to execute the K-means algorithm:

Val MAX_ITERS = 150
def run(K: Int, obs: DblMatrix): Unit = {
  val kmeans = KMeans[Double](K, MAX_ITERS) //1

  val clusters = kmeans |> XTSeries[DblVector](obs) //2
  clusters.foreach( _.center.foreach( show( _ ))) //3
  clusters.map( _.stdDev(XTSeries[DblVector](obs, euclidean))).foreach( show( _ )  )  //4
}

The KMeans class is first initialized with a number of clusters, K, and a maximum number of iterations, MAX_ITERS (line 1). These two parameters are domain and problem specific. The clustering algorithm is executed (line 2) returning a list of clusters. The clusters' centroid information is then displayed (line 3) and the standard deviation is computed for each of the clusters for a given number of clusters, K, and observations, obs (line 4).

Let's load the data from CSV files using the DataSource class (refer to the Data extraction section in Appendix A, Basic Concepts):

final val path = "resources/data/chap4/"
val extractor = YahooFinancials.adjClose :: List[Array[String] =>Double]() // 5
def symbols = DataSource.listSymbols(path)  //6
        
final val START = 80
final val SAMPLES = 50
val normalize=true
val prices = symbols.map(s =>DataSource(s,path,normalize) |> extractor)  //7
prices.find(_.isEmpty) match {  //8
  case Some(noPrice) = { … }
  case None => {
    val values = prices. map(x => x(0))
                        .map(_.drop(START).take(SAMPLES))
    args.map(_.toInt) foreach( run(_, values)) //9
  }
}

As mentioned earlier, the cluster analysis applies to the closing price in the range between the 80^th and 130^th trading day. The extractor is defined to extract the adjusted closing price for a stock whose price information is retrieved from YahooFinancials (line 5). The list of stock symbols is used to extract price information from CSV files located at the path (line 6). For instance, the ticker symbol for General Electric Corp. is GE and the trading data is located in GE.csv.

The 50 daily prices for each stock are extracted by an instance of DataSource (line 7). The run method introduced earlier is invoked either for each stock or as soon as K-means fails through an exit condition in the find method (line 8). The normalized data values.toArray for the specific time window is extracted by the combination of calls to drop and take Scala array methods (line 9).

The first test run is executed with K=3 clusters. The mean (or centroid) vector for each cluster is plotted as follows:

Chart of means of clusters using K-means K=3

The means vectors of the three clusters are quite distinctive. The top and bottom means 1 and 2 in the chart have the respective standard deviation of 0.34 and 0.27 and share a very similar pattern. The difference between the elements of the 1 and 2 cluster mean vectors is almost constant: 0.37. The cluster with a mean vector 3 represents the group of stocks that behave like the stocks in cluster 2 at the beginning of the time period, and behave like the stocks in cluster 1 towards the end of the time period.

This behavior can be easily explained by the fact that the time window or trading period, the 80^th to 130^th trading day, correspond to the shift in the monetary policy of the federal reserve in regard to the quantitative easing program. Here is the list of stocks for each of the clusters whose centroid values are displayed on the chart:

Let's evaluate the impact of the number of clusters K on the characteristics of each cluster.

We repeat the previous test on the 127 stocks and the same time window with the number of clusters varying from 2 to 15.

The mean (or centroid) vector for each cluster is plotted as follows for K = 2:

Chart of means of clusters using K-means K=2

The chart of the results of the K-means algorithms with 2 clusters shows that the mean vector for the cluster labeled 2 is similar to the mean vector labeled 3 on the chart with K = 3 clusters. However, the cluster with the mean vector 1 reflects somewhat the aggregation or summation of the mean vectors for the clusters 1 and 3 in the chart K =3. The aggregation effect explains why the standard deviation for the cluster 1, 0.55, is twice as much as the standard deviation for the cluster 2, 0.28.

The mean (or centroid) vector for each cluster is plotted as follows for K = 5:

Chart of means of clusters using K-means K=5

In this chart, we can assess that the clusters 1 (with the highest mean), 2 (with the lowest mean), and 3 are very similar to the clusters with the same labels in the chart for K =3. The cluster with the mean vector 4 contains stocks whose behaviors are quite similar to those in cluster 3, but in the opposite direction. In other words, the stocks in cluster 3 and 4 reacted in opposite ways following the announcement of the change in the monetary policy.

In the tests with high values of K, the distinction between the different clusters becomes murky, as shown in the following chart for K = 10:

Chart of means of clusters using K-means K=10

The means for clusters 1, 2, and 3 seen in the first chart for the case K = 3 are still visible. It is fair to assume that these are very likely the most reliable clusters. These clusters happened to have a low standard deviation or high density.

Let's define the density of a cluster C_j with a centroid c_j as the inverse of the Euclidean distance between all members of each cluster and its mean (or centroid):

The density of the cluster is plotted against the number of clusters with K = 1 to 13:

Bar chart of the average cluster density for K = 1 to 13

As expected, the average density of each cluster increases as K increases. From this experiment, we can draw the simple conclusion that the density of each cluster does not significantly increase in the test runs for K =5 and beyond. You may observe that the density does not always increase as the number of clusters increases (K = 6 to K = 11). The anomaly can be explained by the following three factors:

There are several methodologies to validate the output of a K-means algorithm from purity to mutual information [4:6]. One effective way to validate the output of a clustering algorithm is to label each cluster and run those clusters through a new batch of labeled observations. For example, if during one of these tests you find that one of the clusters CC contains most of the commodity-related stocks, then you can select another commodity-related stock, SC, which is not part of the first batch, and run the entire clustering algorithm again. If SC is contained in CC, then the clustering has performed as expected. If this is the case, you should run a new set of stocks, some of which are commodity related, and measure the number of true positives, true negatives, false positives, and false negatives. The precision, recall, and F1 measures introduced in the Assessing a model section of Chapter 2, Hello World!, confirms whether the tuning parameters and labels you selected for your cluster are indeed correct.

We reviewed some of the tuning parameters that impact the quality of the results of the K-means clustering. They are as follows:

In some cases, the similarity criterion (that is, Euclidean distance versus cosine value) can have an impact on the cleanness or density of the clusters.

The final and important consideration is the computational complexity of the K-means algorithm. The previous sections of the chapter described some of the performance issues with K-means and possible remedies.

Despite its many benefits, the K-means algorithm does not handle missing data or unobserved features very well. Features that depend on each other indirectly may in fact depend on a common hidden (also known as latent) variable. The expectation-maximization algorithm described in the next section addresses some of these limitations.

Latent variables Z can be visualized as the behavior (or symptoms) of a model (observed) X for which Z are the root causes of the behavior:

Visualization of observed and latent features

The latent values Z follow a Gaussian distribution. For the statisticians among us, the mathematics of a mixture model is described in the following information box.

Note

The mixture model

If {x_i} is a set of observed features associated with latent features {z_k}, the probability for the feature x_i given z_k has a value j:

The probability p is called the base distribution. If we extend to the entire model, θ= {x_i, z_k}, the conditional probability is defined as follows:

The most widely used mixture model is the Gaussian mixture model that represents the base distribution p as a Normal distribution and the conditional probability as a weighted Normal multivariate distribution:

As far as the implementation is concerned, the expectation-maximization algorithm can be broken down into three stages:

A formal, detailed, but short mathematical formulation of the EM algorithm can be found in S. Borman's tutorial [4:9].

Let's implement the three steps (LL, E-step, and M step) in Scala. The internal calculations of the EM algorithm are a bit complex and overwhelming. You may not benefit much from the details of a specific implementation such as computation of the eigenvalues of the covariance matrix of the expectation of the log likelihood. This implementation hides some complexities by using the Apache Commons Math library package [4:10].

First, let's define convenient internal types:

type EM = MultivariateNormalMixtureExpectationMaximization
type EMOutput = List[(Double, DblVector, DblVector)]
import scala.collections.JavaConversions._ //1

The constructor of the MultivariateEM class uses the standard template for machine learning algorithm classes:

Here is an implementation of the constructor of MultivariateEM:

class MultivariateEM[T <% Double](K: Int) extends PipeOperator[XTSeries[Array[T]], EMOutput]

The Apache Commons Math Java implementation of the EM uses Java container classes that need to be explicitly converted to Scala collections. Those conversions are defined in the JavaConversions package (line 1).

The implementation of the EM algorithm in the data transformation |> operator uses the Apache Commons Math MultivariateNormalMixture class for the Gaussian mixture model and the MultivariateNormalMixtureExpectationMaximization class for the EM algorithm:

def |> : PartialFunction[XTSeries[Array[T]], EMOutput] = {
 case xt: XTSeries[Array[T]] if(xt.size>0 && dimension(xt)>0) =>{
   val data: DblMatrix = xt  //2
   val multivariateEM = new EM(data)
   val est = MultivariateNormalMixtureExpectationMaximization
             .estimate(data, K)
   multivariateEM.fit(est) //3

   val newMixture = multivariateEM.getFittedModel //4
   val components = newMixture.getComponents.toList //5
   components.map(p => (p.getKey.toDouble, p.getValue.getMeans, p.getValue.getStandardDeviations)  ))//6
….

Let's look at the main |> method of the MultivariateEM wrapper class. The first step is to convert the time series into a primitive matrix of Double with observations and historical quotes as rows and the stock symbols as columns (line 2).

The initial mixture of Gaussian distributions can be provided by the user or can be extracted from the dataset as an estimate (line 3). The getFittedModel model triggers the M-step (line 4).

The Apache library uses Java primitives that need to be converted to Scala types using the package import scala.collection.JavaConversions. An instance of java.util.List is converted to scala.collection.immutable.List using toList, which invokes the asScalaIterator method of WrapAsScala, one of the base traits of JavaConversions (line 5).

The <Double, MultivariateNormalDistribution> key-value pair, returned by the call to getFittedModel by the Apache math method, is to be converted to a tuple containing the mean and standard deviation for each cluster (line 6).

  Try {
       ..
  } match {
      case Success(results) => …
      case Failure(exception)  => ...
  }

Let's apply the MultivariateEM class to the clustering of the same 127 stocks used in evaluating the K-means algorithm.

As discussed in the paragraph related to the curse of dimensionality, the number of stocks (127) to analyze restricts the number of observations to be used by the EM algorithm. A simple option is to filter out some of the noise of the stocks and apply a basic sampling method. The maximum sampling rate is restricted by the frequencies in the spectrum of noises of different types in the historical price of every stock.

The object of the test is to evaluate the impact of the sampling rate, samplingRate, and the number K of clusters used in the EM algorithm:

val extractor = YahooFinancials.adjClose :: List[Array[String] =>Double]() //1

val period = 8
val samplingRate = 10
val smAv = SimpleMovingAverage[Double](period) //2
val obs = DataSource.listSymbols(path).map(sym => { //3
  val xs = DataSource(sym, path, true) |> extractor //2
  val values : XTSeries[Double] = XTSeries.|>(xs)).head //4
  val filtered = smAv |> values
  filtered.zipWithIndex //5
          .drop(period+1).toArray //6
          .filter( _._2%samplingRate==0)
          .map( _._1)
})

The first step is to extract the historical quotes for all the stocks using the same extractor as in the K-means test case (line 1).

The symbols of the stocks under consideration are extracted from the name of the files in the path directory. The historical data is contained in the CSV file named path/STOCK_NAME.csv (line 3). An implicit conversion is triggered by an assignment of values of the type XTSeries[Double] (line 4). The simple moving average algorithm zeroed out the first period values in the smoothed data, filtered (line 5). Those null values have to be dropped before applying the sampling (line 6).

The first test is to execute the EM algorithm with K=3 clusters and a sampling period of 10 on data smoothed by a simple moving average with a period of 8:

MultivariateEM[Double](K) |> XTSeries[DblVector](obs) foreach (…)

The driver prints the key (line 3), the mean (coordinates of the centroid vector) (line 4), and the standard deviation for each component (cluster).

The sampling of historical prices of the 127 stocks between January 1, 2013 and December 31, 2013 with a frequency of 0.1 hertz produces 24 data points. The following chart displays the mean or centroid of each of the 3 clusters:

Chart of the normalized means per cluster using EM K=3

The mean vectors of clusters 2 and 3 have similar patterns, which may suggest that 2 components or clusters could provide a first insight into the similarity within groups of stocks. The following is a chart of the normalized standard deviation per cluster using EM K = 3:

Chart of the normalized standard deviation per cluster using EM K=3

The distribution of the standard deviation along the mean vector of each cluster can be explained by the fact that the price of stocks from a couple of industries went down in synergy, while others went up as a semihomogenous group following the announcement from the Federal Reserve that the monthly quantity of bonds purchased as part of the quantitative easing program would be reduced in the near future.

Online learning is a powerful strategy for training a clustering model when dealing with very large datasets. This strategy has regained interest from scientists lately. The description of online EM is beyond the scope of this tutorial. However, you may need to know that there are several algorithms available for online EM if you ever have to deal with large datasets: batch EM, stepwise EM, incremental EM, and Monte Carlo EM [4:12].