Added Multi task GP regression model, Laplacian Co Regionalization kernel

Author: mandar2812

dynaml-core/scripts/probModels.scala

@@ -17,7 +17,7 @@ val post = c_model.posterior(350)
 histogram((1 to 2000).map(_ => p.sample()))
 histogram((1 to 2000).map(_ => post.sample()))
 
-val data = new Gaussian(-2.0, 3.25).sample(3000).toStream
+val data = new Gaussian(-2.0, 3.25).sample(2000).toStream
 histogram(data)
 title("Histogram of data")
 
@@ -32,7 +32,7 @@ val iidPrior = IIDRandomVarDistr(prior) _
 scatter(iidPrior(1000).sample())
 hold()
 val likelihood = DataPipe((params: (Double, Double)) =>
-  IIDRandomVarDistr(RandomVariable(new Gaussian(params._1, params._2)))(3000))
+  IIDRandomVarDistr(RandomVariable(new Gaussian(params._1, params._2)))(2000))
 
 val gModel = ProbabilityModel(prior, likelihood)
 

dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/kernels/LaplacianKernel.scala

@@ -45,3 +45,17 @@ class LaplaceCovFunc(private var beta: Double)
   override def gradient(x: Double, y: Double): Map[String, Double] =
     Map("beta" -> 1.0*evaluate(x,y)*math.abs(x-y)/math.pow(state("beta"), 2))
 }
+
+class CoRegLaplaceKernel(bandwidth: Double) extends LocalSVMKernel[Int] {
+
+  override val hyper_parameters: List[String] = List("coRegLB")
+
+  state = Map("coRegLB" -> bandwidth)
+
+  override def gradient(x: Int, y: Int): Map[String, Double] =
+    Map("coRegLB" -> 1.0*evaluate(x,y)*math.abs(x-y)/math.pow(state("coRegLB"), 2))
+
+  override def evaluate(x: Int, y: Int): Double = {
+    math.exp(-1.0*math.abs(x-y)/state("coRegLB"))
+  }
+}

dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/models/gp/MTGPRegressionModel.scala

@@ -0,0 +1,41 @@
+package io.github.mandar2812.dynaml.models.gp
+
+import breeze.linalg.{DenseMatrix, DenseVector}
+import io.github.mandar2812.dynaml.kernels.CovarianceFunction
+import org.apache.log4j.Logger
+
+/**
+  * @author mandar2812 date: 11/8/16.
+  *
+  * Abstract implementation of multi-task gaussian process
+  * as outlined in Lawrence et. al 2012 on arxiv
+  * @tparam I The index set of the GP.
+  */
+class MTGPRegressionModel[I](
+  cov: CovarianceFunction[(I, Int), Double, DenseMatrix[Double]],
+  n: CovarianceFunction[(I, Int), Double, DenseMatrix[Double]],
+  data: Seq[Stream[(I, Double)]],
+  num: Int, numOutputs: Int) extends
+  AbstractGPRegressionModel[
+    Seq[Stream[(I, Double)]],
+    (I, Int)](cov, n, data, num*numOutputs) {
+
+  assert(
+    data.length == numOutputs,
+    "Number of outputs in data should match numOutputs constructor variable"
+  )
+
+  private val logger = Logger.getLogger(this.getClass)
+
+  val noutputs = numOutputs
+
+  /**
+    * Convert from the underlying data structure to
+    * Seq[(I, Y)] where I is the index set of the GP
+    * and Y is the value/label type.
+    **/
+  override def dataAsSeq(data: Seq[Stream[(I, Double)]]): Seq[((I, Int), Double)] =
+  data.zipWithIndex.map((patternSet) =>
+    patternSet._1.map(patternAndLabel => ((patternAndLabel._1, patternSet  ._2), patternAndLabel._2))
+  ).reduceLeft((s1, s2) => s1 ++ s2)
+}

dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/probability/ProbabilityModel.scala

@@ -21,7 +21,9 @@ DistL <: Density[Domain]](
 
   val likelihood: DataPipe[ConditioningSet, RandomVarWithDistr[Domain, DistL]] = c
 
-  var Max_Candidates: Int = 10000
+  var Max_Candidates: Int = 1000
+
+  var Max_Estimations: Int = 100000
 
   override val underlyingDist = new Density[(ConditioningSet, Domain)] {
     override def apply(x: (ConditioningSet, Domain)): Double =
@@ -38,9 +40,9 @@ DistL <: Density[Domain]](
     val sampl = this.prior.sample
     val q = this.prior.underlyingDist
 
-    val M = (1 to Max_Candidates).map(i => {
+    val M = (1 to Max_Estimations).map(i => {
       likelihood(sampl()).underlyingDist(data)
-    }).sum/Max_Candidates.toDouble
+    }).sum/Max_Estimations.toDouble
 
     val postD = new Density[ConditioningSet] {
       override def apply(x: ConditioningSet): Double =