@@ -2,7 +2,7 @@ import breeze.linalg.{DenseMatrix, DenseVector}
22import breeze .stats .distributions .{ContinuousDistr , Gamma }
33import io .github .mandar2812 .dynaml .DynaMLPipe ._
44import io .github .mandar2812 .dynaml .analysis .VectorField
5- import io .github .mandar2812 .dynaml .kernels .{DiracKernel , SEKernel }
5+ import io .github .mandar2812 .dynaml .kernels .{DiracKernel , LaplacianKernel , SEKernel }
66import io .github .mandar2812 .dynaml .modelpipe .GPRegressionPipe
77import io .github .mandar2812 .dynaml .pipes .{DataPipe , StreamDataPipe }
88import io .github .mandar2812 .dynaml .probability .MultGaussianRV
@@ -19,6 +19,7 @@ type Data = Stream[(Features, Features)]
1919implicit val f = VectorField (deltaT)
2020
2121val kernel = new SEKernel (1.5 , 1.5 )
22+ val kernel2 = new LaplacianKernel (2.0 )
2223val noise_kernel = new DiracKernel (1.0 )
2324
2425val data_size = 500
@@ -42,7 +43,7 @@ val prepare_data = {
4243
4344val create_gp_model = GPRegressionPipe (
4445 (d : Data ) => d.toSeq.map(p => (p._1, p._2(0 ))),
45- kernel , noise_kernel
46+ kernel2 , noise_kernel
4647)
4748
4849val model_flow = DataPipe (create_gp_model, scales_flow_stub)
@@ -52,15 +53,21 @@ val workflow = prepare_data > model_flow
5253val (model, scales) = workflow(" data/santafelaser.csv" )
5354
5455val num_hyp = model._hyper_parameters.length
55- val proposal = MultGaussianRV (num_hyp)(DenseVector .zeros[Double ](num_hyp), DenseMatrix .eye[Double ](num_hyp))
56+
57+ val proposal = MultGaussianRV (
58+ num_hyp)(
59+ DenseVector .zeros[Double ](num_hyp),
60+ DenseMatrix .eye[Double ](num_hyp)* 0.001 )
5661
5762val mcmc = HyperParameterMCMC [model.type , ContinuousDistr [Double ]](
58- model, model._hyper_parameters.map(h => (h, new Gamma (1.0 , 1 .0 ))).toMap,
63+ model, model._hyper_parameters.map(h => (h, new Gamma (1.0 , 2 .0 ))).toMap,
5964 proposal)
6065
61- val samples = mcmc.iid(500 ).draw
66+ mcmc.burnIn = 50
67+
68+ val samples = mcmc.iid(1000 ).draw
6269
63- val samples_se = samples.map(h => (h(" bandwidth " ), h(" amplitude " )))
70+ val samples_se = samples.map(h => (h(" beta " ), h(" noiseLevel " )))
6471
6572scatter(samples_se)
6673title(" x,y ~ P(sigma, a | Data)" )
0 commit comments