[R-package][python-package] small change in demo for Vecchia approxim…

…ation

R-package/R/GPModel.R

@@ -997,7 +997,6 @@ gpb.GPModel <- R6::R6Class(
       if (!is.null(vecchia_pred_type)) {
         private$vecchia_pred_type <- vecchia_pred_type
       }
-      vecchia_pred_type_c_str <- private$vecchia_pred_type
       if (!is.null(num_neighbors_pred)) {
         private$num_neighbors_pred <- as.integer(num_neighbors_pred)
       }
@@ -1014,7 +1013,7 @@ gpb.GPModel <- R6::R6Class(
         , gp_coords_pred
         , gp_rand_coef_data_pred
         , X_pred
-        , vecchia_pred_type_c_str
+        , private$vecchia_pred_type
         , private$num_neighbors_pred
         , private$cg_delta_conv_pred
       )
@@ -1057,7 +1056,6 @@ gpb.GPModel <- R6::R6Class(
       if (!is.null(vecchia_pred_type)) {
         private$vecchia_pred_type <- vecchia_pred_type
       }
-      vecchia_pred_type_c_str <- private$vecchia_pred_type
       if (!is.null(num_neighbors_pred)) {
         private$num_neighbors_pred <- as.integer(num_neighbors_pred)
       }
@@ -1336,7 +1334,7 @@ gpb.GPModel <- R6::R6Class(
         , cov_pars
         , X_pred
         , use_saved_data
-        , vecchia_pred_type_c_str
+        , private$vecchia_pred_type
         , private$num_neighbors_pred
         , private$cg_delta_conv_pred
         , fixed_effects

R-package/demo/generalized_linear_Gaussian_process_mixed_effects_models.R

@@ -363,6 +363,13 @@ gp_model <- fitGPModel(gp_coords = coords_train, cov_function = "exponential",
                        vecchia_approx = TRUE, num_neighbors = 30, y = y_train,
                        likelihood = likelihood)
 summary(gp_model)
+# Prediction: setting 'num_neighbors_pred' to a larger value than 'num_neighbors' for training
+#   can lead to better predictions
+pred_vecchia <- predict(gp_model, gp_coords_pred = coords_test, num_neighbors_pred = 100,
+                predict_var = TRUE, predict_response = FALSE)
+ggplot(data = data.frame(s_1=coords_test[,1],s_2=coords_test[,2],b=pred_vecchia$mu),aes(x=s_1,y=s_2,color=b)) +
+  geom_point(size=8, shape=15) + scale_color_viridis(option = "B") + ggtitle("Predicted latent GP mean with Vecchia approxmation")
+
 
 #--------------------Gaussian process model with Wendland covariance function----------------
 gp_model <- fitGPModel(gp_coords = coords_train, cov_function = "wendland", 

examples/python-guide/generalized_linear_Gaussian_process_mixed_effects_models.py

@@ -328,12 +328,14 @@ def simulate_response_variable(lp, rand_eff, likelihood):
 pred_var_plot = pred['var'].reshape((nx, nx))
 CS = axs[1, 0].contourf(coords_test_x1, coords_test_x2, pred_var_plot)
 axs[1, 0].set_title("Predicted latent GP standard deviation")
+plt.show(block=False)
 
 # Predict latent GP at training data locations (=smoothing)
 GP_smooth = gp_model.predict_training_data_random_effects()
 # Compare true and predicted random effects
 plt.scatter(b_train, GP_smooth)
 plt.title("Comparison of true and smoothed GP")
+plt.show(block=False)
 # The above is equivalent to the following
 #GP_smooth2 = gp_model.predict(gp_coords_pred=coords_train)
 #np.sum(np.abs(GP_smooth['GP'] - GP_smooth2['mu']))
@@ -350,6 +352,14 @@ def simulate_response_variable(lp, rand_eff, likelihood):
                        vecchia_approx=True, num_neighbors=30, likelihood=likelihood)
 gp_model.fit(y=y_train)
 gp_model.summary()
+# Prediction: setting 'num_neighbors_pred' to a larger value than 'num_neighbors' for training
+#   can lead to better predictions
+pred_vecchia = gp_model.predict(gp_coords_pred=coords_test, num_neighbors_pred=100,
+                                predict_var=True, predict_response=False)
+pred_vecchia = pred_vecchia['mu'].reshape((nx, nx))
+plt.contourf(coords_test_x1, coords_test_x2, pred_vecchia)
+plt.title("Predicted latent GP mean with Vecchia approxmation")
+plt.show(block=False)
 
 #--------------------Gaussian process model with Wendland covariance function----------------
 gp_model = gpb.GPModel(gp_coords=coords_train, cov_function="wendland",
@@ -375,6 +385,7 @@ def simulate_response_variable(lp, rand_eff, likelihood):
 plt.scatter(b3, GP_smooth['GP_rand_coef_nb_2'], label="2. random coef. GP", alpha=0.5)
 plt.legend()
 plt.title("Comparison of true and smoothed GP")
+plt.show(block=False)
 
 # --------------------Using cluster_ids for independent realizations of GPs----------------
 cluster_ids = np.zeros(ntrain)