$ scala RegressionMonoid , Result: Alpha 3.000 Beta 2.000

RegressionMonoid.scala

import collection.parallel.ParSeq
object RegressionMonoid extends App {

  // Use monoidal addition as opposed to def mean(x:Seq[Double, size:Int) = x.sum/size
  def mean(x:ParSeq[Double], size:Int):Double = {
    def monoid(ab:Seq[Double]) = ab.size match { case 2 => ab.head + ab.last; case _ => ab.head }
    val pairs = x.toList.sliding(2,2).toSeq
    val sumpar = pairs.par.map(monoid)
    if (sumpar.size == 1 ) sumpar.head/size else mean(sumpar, size)
  }

  def product(x:ParSeq[Double], y:ParSeq[Double]) = x.zip(y).map{ab => val(a,b) = ab; a*b}

  /* return alpha, beta
   y = alpha + beta * x
  (Y-ymu)/ysigma = corr * (x-xmu)/xsigma
  since beta = covar = ysigma * corr/xsigma
  So, Y - ymu = (x-xmu) * beta
  Y = x * (covar) + (ymu - xmu*beta)
  alpha = (ymu - xmu*beta)
  */
  def ols(x:Seq[Double], y:Seq[Double]) = {
    val (xnorm, xmu, xsigma) = normalize(x)
    val (ynorm, ymu, ysigma) = normalize(y)
    val xy = product(xnorm, ynorm)
    val corr = mean(xy, xy.size)
    val beta = ysigma * corr/xsigma
    val alpha = ymu - xmu * beta
    (alpha, beta, corr)
  }

  def normalize(x:Seq[Double]) = {
    val xsq = product(x.par,x.par)
    val mu = mean(x.par, x.size)
    val sigma = math.sqrt(mean(xsq, xsq.size) - mu*mu)
    (x.par.map( i=> (i-mu)/sigma), mu, sigma)
  }

  val x = (1.0 to 100000.0 by 0.5).toSeq
  val y = x.map(i => 2*i + 3)
  val (alpha, beta, corr) = ols(x,y)
  printf("Alpha %.3f Beta %.3f Corr %.3f\n", alpha, beta, corr)
}