Santa Fe Laser: Added script for mcmc sampling of GP hyper parameters

Author: mandar2812

scripts/gp_mcmc_santafe.sc

@@ -0,0 +1,68 @@
+import breeze.linalg.{DenseMatrix, DenseVector}
+import breeze.stats.distributions.{ContinuousDistr, Gamma}
+import io.github.mandar2812.dynaml.DynaMLPipe._
+import io.github.mandar2812.dynaml.analysis.VectorField
+import io.github.mandar2812.dynaml.kernels.{DiracKernel, SEKernel}
+import io.github.mandar2812.dynaml.modelpipe.GPRegressionPipe
+import io.github.mandar2812.dynaml.pipes.{DataPipe, StreamDataPipe}
+import io.github.mandar2812.dynaml.probability.MultGaussianRV
+import io.github.mandar2812.dynaml.probability.mcmc.HyperParameterMCMC
+import io.github.mandar2812.dynaml.utils.GaussianScaler
+import com.quantifind.charts.Highcharts._
+
+
+val deltaT = 4
+
+type Features = DenseVector[Double]
+type Data = Stream[(Features, Features)]
+
+implicit val f = VectorField(deltaT)
+
+val kernel = new SEKernel(1.5, 1.5)
+val noise_kernel = new DiracKernel(1.0)
+
+val data_size = 500
+
+val scales_flow_stub = identityPipe[(GaussianScaler, GaussianScaler)]
+
+val prepare_data = {
+  fileToStream >
+    trimLines >
+    extractTrainingFeatures(
+      List(0), Map()
+    ) >
+    DataPipe((lines: Stream[String]) =>
+      lines.zipWithIndex.map(couple => (couple._2.toDouble, couple._1.toDouble))
+    ) >
+    deltaOperation(deltaT, 0) >
+    StreamDataPipe((r: (Features, Double)) => (r._1, DenseVector(r._2))) >
+    gaussianScaling >
+    DataPipe(DataPipe((d: Data) => d.take(data_size)), scales_flow_stub)
+}
+
+val create_gp_model = GPRegressionPipe(
+  (d: Data) => d.toSeq.map(p => (p._1, p._2(0))),
+  kernel, noise_kernel
+)
+
+val model_flow = DataPipe(create_gp_model, scales_flow_stub)
+
+val workflow = prepare_data > model_flow
+
+val (model, scales) = workflow("data/santafelaser.csv")
+
+val num_hyp = model._hyper_parameters.length
+val proposal = MultGaussianRV(num_hyp)(DenseVector.zeros[Double](num_hyp), DenseMatrix.eye[Double](num_hyp))
+
+val mcmc = HyperParameterMCMC[model.type, ContinuousDistr[Double]](
+  model, model._hyper_parameters.map(h => (h, new Gamma(1.0, 1.0))).toMap,
+  proposal)
+
+val samples = mcmc.iid(500).draw
+
+val samples_se = samples.map(h => (h("bandwidth"), h("amplitude")))
+
+scatter(samples_se)
+title("x,y ~ P(sigma, a | Data)")
+xAxis("sigma")
+yAxis("a")