Added variance to examples/julia-first.jl after modifying gp.jl to ha…

…ndle matrices in predict

examples/julia-first.jl

@@ -10,7 +10,7 @@ t = sort(cat(1, 3.8 * rand(57), 5.5 + 4.5 * rand(68)))
 yerr = 0.08 + (0.22-0.08)*rand(length(t))
 y = 0.2*(t-5.0) + sin.(3.0*t + 0.1*(t-5.0).^2) + yerr .* randn(length(t))
 
-true_t = linspace(0, 10, 5000)
+true_t = linspace(0, 10, 1000)
 true_y = 0.2*(true_t-5) + sin.(3*true_t + 0.1*(true_t-5).^2)
 
 PyPlot.clf()
@@ -38,19 +38,24 @@ gp = celerite.Celerite(kernel)
 celerite.compute!(gp, t, yerr)
 celerite.log_likelihood(gp, y)
 
+#mu, variance = celerite.predict_full_ldlt(gp, y, true_t, return_var=true)
 #mu, variance = celerite.predict_full(gp, y, true_t, return_var=true)
-mu  = celerite.predict_full(gp, y, true_t, return_cov=false, return_var=false)
-#sigma = sqrt(variance)
+mu, variance = celerite.predict(gp, y, true_t, return_var=true)
+#mu  = celerite.predict_full(gp, y, true_t, return_cov=false, return_var=false)
+#mu  = celerite.predict!(gp, gp.x, y, true_t)
+sigma = sqrt.(variance)
 
 PyPlot.plot(true_t, true_y, "k", lw=1.5, alpha=0.3,label="Simulated data")
 PyPlot.errorbar(t, y, yerr=yerr, fmt=".k", capsize=0)
 PyPlot.plot(true_t, mu, "g",label="Initial kernel")
-#PyPlot.fill_between(true_t, mu+sigma, mu-sigma, color="g", alpha=0.3)
+PyPlot.fill_between(true_t, mu+sigma, mu-sigma, color="g", alpha=0.3)
 PyPlot.xlabel("x")
 PyPlot.ylabel("y")
 PyPlot.xlim(0, 10)
 PyPlot.ylim(-2.5, 2.5);
+read(STDIN,Char)
 
+PyPlot.clf()
 vector = celerite.get_parameter_vector(gp.kernel)
 mask = ones(Bool, length(vector))
 mask[2] = false  # Don't fit for the first Q
@@ -70,13 +75,16 @@ vector
 celerite.set_parameter_vector!(gp.kernel, vector)
 
 #mu, variance = celerite.predict(gp, y, true_t, return_var=true)
-mu = celerite.predict(gp, y, true_t, return_cov = false, return_var=false)
-#sigma = sqrt(variance)
+#mu, variance = celerite.predict_full_ldlt(gp, y, true_t, return_var=true)
+#mu = celerite.predict(gp, y, true_t, return_cov = false, return_var=false)
+#mu = celerite.predict!(gp, gp.x, y, true_t)
+mu, variance = celerite.predict(gp, y, true_t, return_var=true)
+sigma = sqrt.(variance)
 
 PyPlot.plot(true_t, true_y, "k", lw=1.5, alpha=0.3)
 PyPlot.errorbar(t, y, yerr=yerr, fmt=".k", capsize=0)
 PyPlot.plot(true_t, mu, "r",label="Optimized kernel")
-#PyPlot.fill_between(true_t, mu+sigma, mu-sigma, color="g", alpha=0.3)
+PyPlot.fill_between(true_t, mu+sigma, mu-sigma, color="g", alpha=0.3)
 PyPlot.xlabel("x")
 PyPlot.ylabel("y")
 PyPlot.xlim(0, 10)

src/gp.jl

@@ -786,7 +786,14 @@ function predict_full_ldlt(gp::Celerite, y, t; return_cov=true, return_var=false
 
     KxsT = transpose(Kxs)
     if return_var
-        v = -sum(KxsT .* apply_inverse_ldlt(gp, KxsT), 1)
+#        println("Size of y: ",size(y)," size of t: ",size(t)," size of KxsT: ",size(KxsT))
+        v = zeros(t)
+        for i=1:length(t)
+          v[i] = -sum(KxsT[:,i] .* apply_inverse_ldlt(gp, KxsT[:,i]))
+          if mod(i,100) == 0
+#            println(i," ",v[i])
+          end
+        end
         v = v + get_value(gp.kernel, [0.0])[1]
         return mu, v[1, :]
     end
@@ -808,9 +815,14 @@ function predict(gp::Celerite, y, t; return_cov=true, return_var=false)
 
     KxsT = transpose(Kxs)
     if return_var
-        v = -sum(KxsT .* apply_inverse(gp, KxsT), 1)
+        v=zeros(t)
+        for i=1:length(t)
+#        v = -sum(KxsT .* apply_inverse(gp, KxsT), 1)
+          v[i] = -sum(KxsT[:,i] .* apply_inverse(gp, KxsT[:,i]))
+        end
         v = v + get_value(gp.kernel, [0.0])[1]
-        return mu, v[1, :]
+#        return mu, v[1, :]
+        return mu, v
     end
 
     cov = get_matrix(gp, t)