Main parameters for GPBoost

Gaussian process and random effects model option

Below is a list of parameters for GPModel() objects for modeling Gaussian processes (GPs) and grouped random effects and for specifying how these models are trained.

• Currently supported likelihoods

• Currently supported GP covariance functions including ARD, estimating the smoothness parameter, and space-time models

• Currently supported GP large data approximations such as vecchia and vif approximations

• Optimization parameters for additional optimization options for the params argument of the fit() and set_optim_params() functions including (i) monitoring convergence, (ii) optimization algorithm options, (iii) manually setting initial values for parameters, and (iv) selecting which parameters are estimated. See the the documentation of the Python and R packages for exhaustive lists of all parameters for the params argument.

Model specification parameters

• likelihood : string, (default = gaussian)

  • Likelihood function, i.e., distribution of the response variable conditional on fixed and random effects

  • This is set when defining a GPModel() for both the GPBoost algorithm and (generalized) linear mixed effects and Gaussian process models

  • Currently supported likelihoods, grouped by response type:

    Response type	Support	Likelihoods
    Continuous response	y in (-inf, inf)	gaussian, t, t_fix_df, quantile_regression / asymmetric_laplace, gaussian_heteroscedastic, gaussian_heteroscedastic_fixed_and_random
    Positive continuous response	y in (0, inf)	gamma, gamma_varying_shape, lognormal, gpd, egpd_power, egpd_power_mixture, egpd_beta, egpd_power_beta
    Non-negative continuous / semicontinuous response	y in [0, inf) with point mass at 0	tweedie, tweedie_fixed_p, tweedie_joint, tweedie_varying_dispersion, tweedie_joint_varying_dispersion (and their _fixed_p variants), hurdle_gamma, hurdle_lognormal, hurdle GPD / EGPD likelihoods, zero_censored_power_transformed_normal, zero_censored_shifted_gamma
    Count response	y in {0, 1, 2, ...}	poisson, negative_binomial, negative_binomial_1, zero-inflated count likelihoods
    Binary response	y in {0, 1}	bernoulli_logit, bernoulli_probit
    Proportion / fractional / bounded response	y in [0, 1]	quasi_bernoulli_logit, quasi_bernoulli_probit, binomial_logit, binomial_probit, beta_binomial, beta, zoctn, zero_one_censored_transformed_beta, zero_one_censored_shifted_gamma

    The table lists the main likelihood names. Aliases, heteroscedastic and varying-shape variants, and the _regression variants of the two-part likelihoods are documented in the detailed descriptions below.

    Continuous response: y in (-inf, inf)

    • gaussian : Gaussian likelihood

    • t : t-distribution (e.g., for robust regression). The default approximation is Fisher-Laplace, i.e., the Fisher information is used instead of the observed Hessian in the Laplace approximation of the marginal likelihood.

    • t_fix_df : t-distribution with the degrees-of-freedom (df) held fixed and not estimated

      • The degrees-of-freedom (df) can be set via the likelihood_additional_param parameter. The default is df = 2

    • quantile_regression / asymmetric_laplace : an asymmetric Laplace likelihood for quantile regression. Both names are accepted, as is the alias quantile

      • The quantile must be supplied through likelihood_additional_param and must be strictly between 0 and 1

      • The default approximation is Fisher-Laplace, i.e., the Fisher information is used instead of the observed Hessian in the Laplace approximation of the marginal likelihood.

      • Enable the triangular-kernel-curvature (TKC) approximation by appending _triangular_kernel_curvature or the shorthand _tkc to the likelihood name, for example, quantile_regression_tkc or asymmetric_laplace_triangular_kernel_curvature.

      • experimental: the log-likelihood is not differentiable at the kinks where y_i equals the location parameter, and the quasi-Newton mode finding can therefore stall at a point that is not the exact posterior mode. Appending _ssn_alm, for example quantile_regression_ssn_alm, enables an exact check of the (non-smooth) optimality conditions after the mode finding and, if the check fails, a refinement of the mode with a semismooth Newton method applied to the subproblems of an augmented Lagrangian method (SSN-ALM). This makes the approximate marginal likelihood a well-defined function of the parameters, i.e., independent of the path taken by the optimizer and of the number of threads, at the price of additional linear solves during mode finding. This matters for standard errors and for likelihood-based model comparison, but it does not necessarily improve predictive accuracy; increasing max_num_restarts_lbfgs can be a cheaper way of avoiding poor optima. Appending _admm_ssn_alm instead warm starts the refinement with an ADMM phase, which typically halves its cost when matrix_inversion_method = "cholesky" and is not recommended for iterative methods.

    • gaussian_heteroscedastic : Gaussian likelihood where the mean is related to fixed and random effects and the log-error variance is related to fixed effects only (linear predictor or GPBoost algorithm). The estimated coefficients of the log-variance model are returned alongside the mean-model coefficients (with the suffix ‘_scale’). Fisher-Laplace is the default and currently the only implemented approximation.

    • gaussian_heteroscedastic_fixed_and_random : Gaussian likelihood where both the mean and the variance are related to fixed and random effects. This is currently only implemented for GPs with a vecchia approximation. Fisher-Laplace is the default and currently the only implemented approximation.

    Positive continuous response: y in (0, inf)

    • gamma : Gamma likelihood with a log link function

    • gamma_varying_shape : As gamma, but the gamma shape varies across observations and is modeled by an additional fixed-effects-only predictor: log(shape) = F_s(X) (F_s(X) = linear predictor or the GPBoost algorithm), while the log mean log(mu) = F(X) + Zb is related to both fixed and random effects. Note that the shape also governs the dispersion, var(y) = mu^2 / shape. The estimated coefficients of the log-shape model are returned alongside the mean-model coefficients (with the suffix ‘_shape’). See also the corresponding two-part likelihoods hurdle_gamma_varying_shape and hurdle_regression_gamma_varying_shape below

    • lognormal : Log-normal likelihood with a log link function

    • gpd : Generalized Pareto likelihood. The log scale parameter equals the latent predictor ‘eta’ (sum of fixed and random effects), ‘sigma = exp(eta)’, and the estimated auxiliary parameter is ‘shape’ with the regular domain ‘shape > -0.5’

    • egpd_power : Naveau extended generalized Pareto likelihood with carrier ‘G(u) = u^kappa’ and auxiliary parameters ‘shape’ and ‘kappa’

    • egpd_power_mixture : Naveau power-mixture carrier with ordered exponents ‘kappa2 = kappa1 + delta_kappa’ and auxiliary parameters ‘shape’, ‘kappa1’, ‘delta_kappa’, and ‘p’. Both exponent parameters are positive and ‘0 < p < 1’

    • egpd_beta : Naveau beta-carrier extended generalized Pareto likelihood with auxiliary parameters ‘shape’ and ‘delta’

    • egpd_power_beta : Naveau power-beta carrier with auxiliary parameters ‘shape’, ‘delta’, and ‘kappa’

    Non-negative continuous / semicontinuous response: y in [0, inf) with point mass at 0

    • tweedie : Compound Poisson–Gamma Tweedie likelihood with a log link, where the latent predictor is ‘eta’, the mean is ‘mu = exp(eta)’, and ‘Var(y | eta) = phi * mu^p’, with ‘1.01 < p < 1.99’. Both dispersion ‘phi’ and power ‘p’ are estimated

    • tweedie_fixed_p : The same Tweedie likelihood with ‘p’ fixed through likelihood_additional_param and only ‘phi’ estimated. The fixed power is mandatory and must satisfy ‘1.01 < p < 1.99’. Fits at different fixed powers include the complete density and can therefore be compared by marginal log-likelihood for power profiling

    • tweedie_joint, tweedie_joint_fixed_p : The same Tweedie model, but the observed number of events (e.g., claims) ‘N’, given in the first column of additional_likelihood_data, is used jointly with the aggregate response ‘y’. The likelihood is the joint density of ‘(y, N)’ of the compound Poisson–Gamma representation, ‘N ~ Poisson(mu^(2-p) / (phi * (2-p)))’ and ‘y | N = n ~ Gamma(n * (2-p) / (p-1), scale = phi * (p-1) * mu^(p-1))’, see Jorgensen and de Souza (1994, Scandinavian Actuarial Journal) “Fitting Tweedie’s compound Poisson model to insurance claims data”. Given ‘phi’ and ‘p’, the mean model is the same as for tweedie; ‘N’ adds information for estimating ‘phi’ and ‘p’. ‘N’ must be a non-negative integer that is 0 if and only if ‘y’ is 0. Predictions are the same as for tweedie, and the test_neg_log_likelihood metric uses the marginal density of ‘y’ since ‘N’ is not known for new data

    • tweedie_varying_dispersion, tweedie_varying_dispersion_fixed_p : As tweedie and tweedie_fixed_p, but the dispersion ‘phi’ varies across observations: ‘log(phi) = F_d(X)’ is related to fixed effects only (linear predictor or GPBoost algorithm), while ‘log(mu) = F(X) + Zb’ is related to both fixed and random effects. The power ‘p’ is then the only (auxiliary) parameter. The estimated coefficients of the log-dispersion model are returned alongside the mean-model coefficients (with the suffix ‘_dispersion’)

    • tweedie_joint_varying_dispersion, tweedie_joint_varying_dispersion_fixed_p : The joint ‘(y, N)’ likelihood of tweedie_joint with the varying dispersion of tweedie_varying_dispersion

    • hurdle_<base> : Two-part likelihoods for non-negative response variables with an excess probability ‘p0’ of exact zeros. They combine a point mass ‘p0’ at zero with a base distribution with support ‘y > 0’ for the remaining probability mass ‘1 - p0’. The fixed effects ‘F(X)’ and random effects ‘Zb’ enter only through the base component: ‘exp(F(X) + Zb)’ is its mean or scale parameter, not the unconditional response mean, which is ‘E(y) = (1 - p0) * base_mean’. The structural-zero probability ‘p0’ is estimated jointly with the base auxiliary parameters. Currently supported bases: hurdle_gamma, hurdle_lognormal, and the extreme-value bases hurdle_gpd, hurdle_egpd_power, hurdle_egpd_power_mixture, hurdle_egpd_beta and hurdle_egpd_power_beta (for these ‘exp(F(X) + Zb)’ is the GPD/EGPD scale parameter and response moments exist only for small enough shape). For all these bases, the prefix zero_inflated_ is accepted as an alias, for example zero_inflated_gamma maps to hurdle_gamma. See zero_inflated_<base> under count responses for the corresponding count-data family

    • hurdle_regression_<base> : As hurdle_<base>, but the structural-zero probability is modeled as a logistic regression on the covariates ‘X’ instead of as a constant (e.g., hurdle_regression_gamma, hurdle_regression_lognormal, hurdle_regression_gpd, hurdle_regression_egpd_power, hurdle_regression_egpd_power_mixture, hurdle_regression_egpd_beta, hurdle_regression_egpd_power_beta). The structural-zero probability is then ‘pi_i = 1 / (1 + exp(-x_i^T alpha))’, modeled through a second fixed-effects-only predictor that reuses the same design matrix ‘X’ as the response model; the response predictor ‘eta’ carries the random effects while the zero predictor does not. The estimated zero-model coefficients ‘alpha’ are returned alongside the response-model coefficients (with the suffix ‘_zero’)

    • hurdle_gamma_varying_shape, hurdle_regression_gamma_varying_shape : As hurdle_gamma and hurdle_regression_gamma, but with the observation-specific gamma shape of gamma_varying_shape. For hurdle_regression_gamma_varying_shape there are three predictors: the response mean, the structural-zero logit (coefficients with the suffix ‘_zero’), and log(shape)

    • zero_censored_power_transformed_normal : Likelihood of a censored and power-transformed normal variable for modeling data with a point mass at 0 and a continuous distribution for y > 0. The model used is Y = max(0,X)^lambda, X ~ N(mu, sigma^2), where mu = F(X) + Zb, and sigma and lambda are (auxiliary) parameters that are estimated. For more details on this model, see Sigrist et al. (2012, AOAS) “A dynamic nonstationary spatio-temporal model for short term prediction of precipitation”

    • zero_censored_power_transformed_normal_heteroscedastic : As zero_censored_power_transformed_normal, but the standard deviation sigma of the latent normal variable varies across observations: log(sigma) = F_2(X) is related to fixed effects only (linear predictor or GPBoost algorithm), while mu = F(X) + Zb is related to both fixed and random effects. lambda is then the only (auxiliary) parameter that is estimated. The estimated coefficients of the log-sigma model are returned alongside the mean-model coefficients (with the suffix ‘_scale’)

    • zero_censored_shifted_gamma : Zero-censored shifted gamma likelihood for modeling data with a point mass at 0 and a continuous distribution for y > 0. The model used is Y = max(Z - xi, 0), where Z follows a gamma distribution with mean mu = exp(F(X) + Zb) and shape k. The shape k and shift xi are (auxiliary) parameters that are estimated. This is the version of zero_one_censored_shifted_gamma (Sigrist and Stahel, 2011) without the upper censoring at 1

    • zero_censored_shifted_gamma_varying_shape : As zero_censored_shifted_gamma, but the shape k varies across observations: log(k) = F_s(X) is related to fixed effects only (linear predictor or GPBoost algorithm), while log(mu) = F(X) + Zb is related to both fixed and random effects. The shift xi is then the only (auxiliary) parameter that is estimated. The estimated coefficients of the log-shape model are returned alongside the mean-model coefficients (with the suffix ‘_shape’)

    Count response: y in {0, 1, 2, …}

    • poisson : Poisson likelihood with log link function

    • negative_binomial : Negative binomial likelihood with a log link function (aka nbinom2, negative_binomial_2). The variance is mu * (mu + r) / r, mu = mean, r = shape, with this parametrization

    • negative_binomial_1 : Negative binomial 1 (aka nbinom1) likelihood with a log link function. The variance is mu * (1 + phi), mu = mean, phi = dispersion, with this parametrization

    • zero_inflated_<base> : Two-part count likelihoods that combine a point mass ‘p0’ at zero with a count base distribution for the remaining probability mass ‘1 - p0’ (the base can itself generate additional zeros). As for the hurdle likelihoods, ‘exp(F(X) + Zb)’ is the mean of the (non-structural) base component - not the unconditional response mean - so that ‘E(y) = (1 - p0) * base_mean’, and the structural-zero probability ‘p0’ is estimated jointly with the base auxiliary parameters. Currently supported bases: zero_inflated_poisson, zero_inflated_negative_binomial (auxiliary parameter: shape; aliases zero_inflated_nbinom2, zero_inflated_negative_binomial_2) and zero_inflated_negative_binomial_1 (auxiliary parameter: dispersion; alias zero_inflated_nbinom1)

    • zero_inflated_regression_<base> : As zero_inflated_<base>, but with the logistic structural-zero model described under hurdle_regression_<base> above: zero_inflated_regression_poisson, zero_inflated_regression_negative_binomial, and zero_inflated_regression_negative_binomial_1

    Binary response: y in {0, 1}

    • bernoulli_logit : Bernoulli likelihood with a logit link function for binary classification. Aliases: binary, binary_logit

    • bernoulli_probit : Bernoulli likelihood with a probit link function for binary classification. Aliases: binary_probit

    Proportion / fractional / bounded response: y in [0, 1]

    • quasi_bernoulli_logit : quasi-Bernoulli likelihood with a logit link function for y in [0,1]. Aliases: quasi_binary, quasi_binary_logit

    • quasi_bernoulli_probit : quasi-Bernoulli likelihood with a probit link function for y in [0,1]. Aliases: quasi_binary_probit

    • binomial_logit : Binomial likelihood with a logit link function. The response variable y needs to contain proportions of successes / trials, and the weights parameter needs to contain the numbers of trials. Aliases: binomial

    • binomial_probit : Binomial likelihood with a probit link function. The response variable y needs to contain proportions of successes / trials, and the weights parameter needs to contain the numbers of trials

    • beta_binomial : Beta-binomial likelihood with a logit link function. The response variable y needs to contain proportions of successes / trials, and the weights parameter needs to contain the numbers of trials. Aliases: betabinomial, beta-binomial

    • beta : Beta likelihood with a logit link function (parametrization of Ferrari and Cribari-Neto, 2004). The response must lie strictly in (0,1)

    • zoctn : Zero-one censored transformed normal likelihood for modeling data in [0,1] with point masses at 0 and 1 and a continuous distribution on (0,1). The model used is T ~ N(mu, sigma^2), W = max(min(T,1),0), and Y = g(W), where g(x) = expit(a + b * logit(x)) for x in (0,1), mu = F(X) + Z_RE u, u denotes the random effects, Z_RE is their design matrix, and sigma, a, and b are (auxiliary) parameters that are estimated. For more details on this model, see Qiang and Sigrist (2026)

    • zero_one_censored_transformed_beta : Zero-one censored transformed beta likelihood for modeling data in [0,1] with point masses at 0 and 1 and a continuous distribution on (0,1). If T follows a beta distribution with mean mu = expit(F(X) + Zb) and precision phi, the observed response is obtained by applying the linear transformation Y = (1 + 2u) * T - u and censoring the result to [0,1]. The precision phi and shift u are (auxiliary) parameters that are estimated. For more details on this model, see Kosmidis and Zeileis (2025)

    • zero_one_censored_shifted_gamma : Zero-one censored shifted gamma likelihood for modeling data in [0,1] with point masses at 0 and 1 and a continuous distribution on (0,1). The model used is Y = min(max(Z - xi, 0), 1), where Z follows a gamma distribution with mean mu = exp(F(X) + Zb) and shape k. The shape k and shift xi are (auxiliary) parameters that are estimated. For more details on this model, see Sigrist and Stahel (2011)

    Notes

    • The first lines in the likelihoods source file contain additional comments on the specific parametrizations used

    • Other likelihoods can be implemented upon request

• group_data : two dimensional array / matrix of doubles or strings, optional (default = None)

  • Labels of group levels for grouped random effects

• group_rand_coef_data : two dimensional array / matrix of doubles or None, optional (default = None)

  • Covariate data for grouped random coefficients

• ind_effect_group_rand_coef : integer vector / array of integers or None, optional (default = None)

  • Indices that relate every random coefficients to a “base” intercept grouped random effect. Counting starts at 1.

• gp_coords : two dimensional array / matrix of doubles or None, optional (default = None)

  • Coordinates (input features) for Gaussian process

• gp_rand_coef_data : two dimensional array / matrix of doubles or None, optional (default = None)

  • Covariate data for Gaussian process random coefficients

• cov_function : string, (default = exponential)

  • Covariance function for the Gaussian process. Available options:

    • matern : Matern covariance function with the smoothness specified by the cov_fct_shape parameter (using the parametrization of Rasmussen and Williams, 2006)

    • matern_estimate_shape : same as matern but the smoothness parameter is also estimated

    • matern_space_time : Spatio-temporal Matern covariance function with different range parameters for space and time

      • Note that the first column in gp_coords must correspond to the time dimension

    • space_time_gneiting : Spatio-temporal covariance function given in Eq. (16) of Gneiting (2002)

      • Note that the first column in gp_coords must correspond to the time dimension

      • This covariance has seven parameters (in the following order: sigma2, a, c, alpha, nu, beta, delta) which are all estimated by default. You can disable the estimation of some of these parameter using the estimate_cov_par_index argument of the params argument in either the fit function of a gp_model object or the set_optim_params function prior to estimation

    • matern_ard: Anisotropic Matern covariance function with Automatic Relevance Determination (ARD), i.e., with a different range parameter for every coordinate of gp_coords

    • matern_ard_estimate_shape : same as matern_ard but the smoothness parameter is also estimated

    • exponential : Exponential covariance function (using the parametrization of Diggle and Ribeiro, 2007)

    • gaussian : Gaussian, aka squared exponential, covariance function (using the parametrization of Diggle and Ribeiro, 2007)

    • gaussian_ard: Anisotropic Gaussian, aka squared exponential, covariance function with Automatic Relevance Determination (ARD), i.e., with a different range parameter for every coordinate of gp_coords

    • powered_exponential : Powered exponential covariance function with the exponent specified by cov_fct_shape parameter (using the parametrization of Diggle and Ribeiro, 2007)

    • wendland : Compactly supported Wendland covariance function (using the parametrization of Bevilacqua et al., 2019, AOS)

    • linear: Linear covariance function. This corresponds to a Bayesian linear regression model with a Gaussian prior on the coefficients with a constant variance diagonal prior covariance, and the prior variance is estimated using empirical Bayes.

    • hurst: Hurst covariance function cov(s, s’) = (sigma2 / 2) * ( ||s||^(2H) + ||s’||^(2H) - ||s - s’||^(2H) ), 0 < H < 1. For H = 0.5, this corresponds to Brownian motion (-> see the estimate_cov_par_index argument). This is the covariance for cov_fct_order = 1 (default). For cov_fct_order = 2, the second-order Hurst covariance cov(s, s’) = sigma2 / (2 (2H - 1)) * ( ||s - s’||^(2H) - ||s||^(2H) - ||s’||^(2H) + 2H (s^T s’) (||s||^(2H-2) + ||s’||^(2H-2)) ), 1 < H < 2, is used. For H = 1.5, this corresponds to an integrated Wiener process. See cov_fct_order for more details

    • hurst_ard: Hurst covariance function with with Automatic Relevance Determination (ARD), i.e., with a different range parameter for every coordinate of gp_coords except for the first coordinate which has a range parameter of 1 due to identifiability with the marginal variance: cov(s, s’) = (sigma2 / 2) * ( (s_1^2 + sum_{k=2}^d (s_k / l_k)^2)^H + (s’_1^2 + sum_{k=2}^d (s’_k / l_k)^2)^H - ((s_1 - s’_1)^2 + sum_{k=2}^d ((s_k - s’_k) / l_k)^2)^H ). For cov_fct_order = 2, the second-order Hurst covariance (see hurst) is applied to the scaled coordinates (s_1, s_2/l_2, …, s_d/l_d)

    • ar1_mf_<base>: Two-level autoregressive multifidelity covariance constructed from a supported base covariance function <base>. For example, use ar1_mf_matern, ar1_mf_matern_ard, or ar1_mf_matern_estimate_shape.

      • The last column of gp_coords is the fidelity indicator and must equal 0 for low-fidelity observations and 1 for high-fidelity observations. All preceding columns are input coordinates for the GPs.

      • The model is

        \[f_H(x) = \rho f_L(x) + \delta(x),\]

        where \(f_L\) and \(\delta\) are independent Gaussian processes with the same covariance-function type but separate parameter vectors.

      • The covariance parameters are ordered as [low-fidelity base parameters, discrepancy base parameters, rho]. The two base-parameter blocks follow the ordinary ordering of <base>. rho is unrestricted and can be negative.

      • All supported base covariance functions except wendland can be used. Correlation tapering and Gaussian-process random coefficients are currently not supported for ar1_mf_<base>.

      • By default, marginal means are fidelity-specific (fidelity_specific_mean = true). For linear regression, a supplied design matrix \(X\) is internally replaced by \([X I(s=0), X I(s=1)]\), yielding independent low- and high-fidelity coefficient vectors. Coefficients are ordered as the low-fidelity block followed by the high-fidelity block. For the GPBoost algorithm, the fidelity indicator is automatically appended to the boosting features in training and prediction, allowing the tree mean to differ by fidelity. Prediction coordinates including the fidelity indicator must therefore be supplied. Set fidelity_specific_mean = false to use one shared marginal mean.

• cov_fct_shape : double, (default = 1.5)

  • Shape parameter of the covariance function (e.g., smoothness parameter for Matern and Wendland covariance). This parameter is irrelevant for some covariance functions such as the exponential or Gaussian.

• cov_fct_order : integer, (default = 1)

  • Order m of the hurst and hurst_ard covariance functions (also when used as base covariance in ar1_mf_hurst and ar1_mf_hurst_ard). Currently, the orders 1 and 2 are supported. For all other covariance functions, this must be 1.

  • The Hurst exponent H of the order m satisfies m - 1 < H < m, and H = m - 0.5 is used as initial value. The covariance parameters are the same for all orders (sigma2, H, and, for hurst_ard, the ranges).

  • Order 1 is a fractional Brownian motion / field, and order 2 is a second-order Hurst process / field (in one dimension, an m-th order fractional Brownian motion, Perrin et al., 2001). Both are anchored at the origin: b(0) = 0 for order 1, and b(0) = 0 and grad b(0) = 0 for order 2. The origin of the coordinates is thus a part of the model, and the coordinates are not centered internally. Shift the coordinates if another anchor is more meaningful.

  • The order is a structural choice and is not estimated. For every order, H is estimated in (m - 1, m), and different orders can be compared, e.g., using the marginal likelihood or cross-validation.

  • H determines the behavior of the process at small scales. If the error variance (or, for non-Gaussian likelihoods, the noise of the observations) is large compared to the variation of the process between neighboring points, H is only weakly identified, and its estimate can be close to the boundaries m - 1 or m. Fixing H (e.g., to m - 0.5, see below) can then be preferable.

  • Continuous-time RW1 and RW2 priors: for a one-dimensional time coordinate t >= 0 whose origin t = 0 is the desired anchor (shift the time points if necessary):

    • Continuous-time RW1: cov_function = "hurst", cov_fct_order = 1, and H fixed to 0.5. This is a Brownian motion with b(0) = 0 and b’(t) = sqrt(q) W’(t), where W’ denotes white noise and q = sigma2.

    • Continuous-time RW2: cov_function = "hurst", cov_fct_order = 2, and H fixed to 1.5. This is an integrated Wiener process with b(0) = b’(0) = 0 and b’’(t) = sqrt(q) W’(t), where q = 3 * sigma2.

    • H is fixed using the init_cov_pars and estimate_cov_par_index parameters. For instance, for a Gaussian likelihood, the covariance parameters are (error variance, sigma2, H), and the following fits a model with a continuous-time RW2 prior:

      X <- cbind(1, time) # unpenalized intercept and linear trend
      gp_model <- fitGPModel(gp_coords = time, cov_function = "hurst", cov_fct_order = 2,
                             likelihood = "gaussian", y = y, X = X,
                             params = list(init_cov_pars = c(0.1, 1, 1.5),
                                           estimate_cov_par_index = c(1, 1, 0)))
      

      X = np.column_stack((np.ones(len(time)), time)) # unpenalized intercept and linear trend
      gp_model = gpb.GPModel(gp_coords=time, cov_function="hurst", cov_fct_order=2, likelihood="gaussian")
      gp_model.fit(y=y, X=X, params={"init_cov_pars": np.array([0.1, 1., 1.5]),
                                     "estimate_cov_par_index": np.array([1, 1, 0])})
      

    • The anchored covariance fixes the polynomial null space of the corresponding intrinsic model (constants for RW1, constants and linear functions for RW2) at the origin. If these components should not be penalized, include them as fixed effects: an intercept for RW1, and an intercept and a linear time trend for RW2.

    • The order 2 model with H = 1.5 is the continuous-time RW2 (integrated Wiener) prior evaluated at the time points, and not the conventional discrete intrinsic RW2 prior on a lattice with independent second differences. For equally spaced time points with spacing h, the second differences of the continuous-time RW2 have variance 2/3 q h^3, and adjacent second differences have correlation 1/4, whereas they are independent for the discrete RW2. The two priors have the same dominant cubic generalized covariance but differ at the scale of the lattice. For multidimensional coordinates, the term “higher-order Hurst field” is less ambiguous than RW1 / RW2.

• gp_approx : string, (default = none)

  • Specifies the use of a large data approximation for Gaussian processes. Available options:

    • none : No approximation

    • vecchia : Vecchia approximation; see Sigrist (2022, JMLR) for more details

      • For matern_space_time and the anisotropic ARD covariance functions (matern_ard, gaussian_ard, matern_ard_estimate_shape), neighbors are selected according to the largest absolute correlations, i.e., using distances in coordinates scaled by the range parameters, and they are redetermined during parameter estimation.

      • For space_time_gneiting and ar1_mf_<base>, neighbors are also selected according to the largest absolute correlations by default. Use gp_approx = vecchia_euclidean for Euclidean-distance selection.

    • full_scale_vecchia : Vecchia-inducing points full-scale (VIF) approximation; see Gyger, Furrer, and Sigrist (2025) for more details

    • tapering : The covariance function is multiplied by a compactly supported Wendland correlation function

    • fitc: Fully Independent Training Conditional approximation aka modified predictive process approximation; see Gyger, Furrer, and Sigrist (2024) for more details

    • full_scale_tapering: Full-scale approximation combining an inducing point / predictive process approximation with tapering on the residual process; see Gyger, Furrer, and Sigrist (2024) for more details

• cluster_ids : one dimensional array (vector) with integer data or Null, (default = Null)

  • IDs / labels indicating independent realizations of random effects / Gaussian processes (same values = same process realization)

• weights : one dimensional array (vector) with numeric data or Null, (default = Null)

  • Sample weights. For a Gaussian likelihood, the error variance (“nugget”) for observation i is divided by weights[i]. For non-Gaussian likelihoods, the conditional log-likelihood contribution of observation i is multiplied by weights[i]. Consequently, weights affect the estimation of both random and fixed effects. Note that a Gaussian likelihood is calculated via the Laplace approximation when gp_approx = "vecchia_latent", when likelihood = "gaussian_latent", and when grouped random effects are combined with a Vecchia-approximated Gaussian process. In these cases, the weights act as for a non-Gaussian likelihood, i.e., the Gaussian log-likelihood contribution of observation i is multiplied by weights[i] instead of the error variance being divided by it.

• additional_likelihood_data : one or two dimensional array (vector or matrix) with numeric data or Null, (default = Null)

  • Observation-level data that some likelihoods require in addition to the response variable y (one row per data point). Currently, this is only used by the joint Tweedie likelihoods (tweedie_joint and its variants), for which its first column contains the observed number of events (e.g., claims) ‘N’. ‘N’ must be a non-negative integer that is 0 if and only if ‘y’ is 0

• cov_fct_taper_range : double, (default = 1.)

  • Range parameter of the Wendland covariance function and Wendland correlation taper function. We follow the notation of Bevilacqua et al. (2019, AOS)

• cov_fct_taper_shape : double, (default = 1.)

  • Shape parameter of the Wendland covariance function and Wendland correlation taper function. We follow the notation of Bevilacqua et al. (2019, AOS)

• num_neighbors : integer

  • Number of neighbors for the Vecchia approximation

  • Internal default values if None:

    • 20 for gp_approx = vecchia

    • 30 for gp_approx = full_scale_vecchia

• vecchia_ordering : string, (default = random)

  • Ordering used in the Vecchia approximation. Available options:

    • none: the default ordering in the data is used

    • random: a random ordering

    • time: ordering accorrding to time (only for space-time models)

    • time_random_space: ordering according to time and randomly for all spatial points with the same time points (only for space-time models)

• vecchia_pred_type : string, (default = Null)

  • Type of Vecchia approximation used for making predictions

  • Default value if vecchia_pred_type = Null : order_obs_first_cond_obs_only

  • Available options:

    • order_obs_first_cond_obs_only : observed data is ordered first and the neighbors are only observed points

    • order_obs_first_cond_all : observed data is ordered first and the neighbors are selected among all points (observed + predicted)

    • latent_order_obs_first_cond_obs_only : Vecchia approximation for the latent process and observed data is ordered first and neighbors are only observed points

    • latent_order_obs_first_cond_all : Vecchia approximation for the latent process and observed data is ordered first and neighbors are selected among all points

    • order_pred_first : predicted data is ordered first for making predictions. This option is only available for Gaussian likelihoods

• num_neighbors_pred : integer, (default = Null)

  • Number of neighbors for the Vecchia approximation for making predictions.

  • Default value if num_neighbors_pred = Null: num_neighbors_pred = 2 * num_neighbors

• num_ind_points : integer

  • Number of inducing points / knots for FITC, full_scale_tapering, and VIF approximations.

  • Internal default values if None:

    • 500 for gp_approx = FITC and gp_approx = full_scale_tapering

    • 200 for gp_approx = full_scale_vecchia

• matrix_inversion_method : string, (default = cholesky)

  • Method used for inverting covariance matrices. Available options:

    • cholesky : Cholesky factorization

    • iterative : iterative methods. A combination of the conjugate gradient, Lanczos algorithm, and other methods.

      This is currently only supported for the following cases:

      • grouped random effects with more than one level

      • likelihood != gaussian and gp_approx == vecchia (non-Gaussian likelihoods with a Vecchia-Laplace approximation)

      • likelihood != gaussian and gp_approx == full_scale_vecchia (non-Gaussian likelihoods with a VIF approximation)

      • likelihood == gaussian and gp_approx == full_scale_tapering (Gaussian likelihood with a full-scale tapering approximation)

• seed : integer, (default = 0)

  • The seed used for model creation (e.g., random ordering in Vecchia approximation)

Optimization parameters

The following list shows some options for the parameter optimization GPModel objects (containing Gaussian process and/or grouped random effects models). These parameters are passed to the params argument of either the fit() function of a GPModel object or to the set_optim_params() function prior to running the GPBoost algorithm. See the the documentation of the Python and R packages for exhaustive lists of all parameters for the params argument.

• trace : bool, optional (default = False)

  • If True, information on the progress of the parameter optimization is printed.

• init_cov_pars : numeric vector / array of doubles, optional (default = Null)

  • Initial values for covariance parameters of Gaussian process and random effects (can be Null). The order it the same as the order of the parameters in the summary function: first is the error variance (only for gaussian likelihood), next follow the variances of the grouped random effects (if there are any, in the order provided in ‘group_data’), and then follow the marginal variance and the range of the Gaussian process. If there are multiple Gaussian processes, then the variances and ranges follow alternatingly. If ‘init_cov_pars = Null’, an internatl choice is used that depends on the likelihood and the random effects type and covariance function. If you select the option ‘trace = true’ in the ‘params’ argument, you will see the first initial covariance parameters in iteration 0.

• init_coef : numeric vector / array of doubles, optional (default = Null)

  • Initial values for the regression coefficients (if there are any, can be Null)

• init_aux_pars : numeric vector / array of doubles, optional (default = Null)

  • Initial values for additional parameters for non-Gaussian likelihoods (e.g., shape parameter of a gamma or negative binomial likelihood) (can be None).

• estimate_cov_par_index : numeric vector / array of integers or NULL, optional (default = -1)

  • This allows for disabling the estimation of some (or all) covariance parameters. If estimate_cov_par_index = -1, all covariance parameters are estimated. If estimate_cov_par_index != -1, this should be a vector with length equal to the number of covariance parameters, and estimate_cov_par_index[i] should be of bool type indicating whether parameter number i is estimated or not. For instance, estimate_cov_par_index = [1,1,0] means that the first two covariance parameters are estimated and the last one not.

  • Parameters that are not estimated are kept at their initial values (see init_cov_pars). They are treated as known constants when calculating the standard errors of the estimated parameters, and their own standard errors are NaN.

• estimate_aux_pars: bool, (default = True)

  • If True, any additional parameters for non-Gaussian likelihoods are also estimated (e.g., shape parameter of a gamma or negative binomial likelihood)

• optimizer_cov : string, optional (default = lbfgs for linear mixed effects models and gradient_descent for the GPBoost algorithm)

  • Optimizer used for estimating covariance parameters

  • Options: “lbfgs”, “gradient_descent”, “fisher_scoring”, “newton” ,”nelder_mead”

  • If there are additional auxiliary parameters for non-Gaussian likelihoods, ‘optimizer_cov’ is also used for those

• optimizer_coef : string, optional (default = wls for Gaussian data and lbfgs for other likelihoods)

  • Optimizer used for estimating linear regression coefficients, if there are any (for the GPBoost algorithm there are usually none)

  • Options: gradient_descent, lbfgs, wls, nelder_mead. Gradient descent steps are done simultaneously with gradient descent steps for the covariance paramters. wls refers to doing coordinate descent for the regression coefficients using weighted least squares

  • If optimizer_cov is set to nelder_mead or lbfgs, optimizer_coef is automatically also set to the same value

• maxit : integer, optional (default = 1000)

  • Maximal number of iterations for optimization algorithm

• delta_rel_conv : double, optional (default = 1e-6 except for nelder_mead for which the default is 1e-8)

  • Convergence tolerance. The algorithm stops if the relative change in eiher the (approximate) log-likelihood or the parameters is below this value.

  • If delta_rel_conv = -999, internal default values are used (= 1e-6 except for nelder_mead for which the default is 1e-8)

Options for the GPBoost algorithm

Metrics for parameter tuning

It is important that tuning parameters (= hyperparameters) for the tree-boosting part are chosen appropriately. There are no universal good “default” values for different data sets. See below for a list of important tuning parameters. Selecting tuning parameters can be done conveniently via the gpb.grid.search.tune.parameters function in the Python and R packages.

The metric parameter (e.g., for the gpb.train, gpboost, and gpb.grid.search.tune.parameters functions in R and Python) specifies how prediction accuracy is measured on validation data.

• For the GPBoost algorithm, i.e., if there is a gp_model, test_neg_log_likelihood is the default metric.

• Other supported metrics include: mse, rmse, mae, crps_gaussian, binary_logloss, binary_error, and auc.

• If another metric besides test_neg_log_likelihood is used for the GPBoost algorithm, it is calculated as follows. First, the predictive mean of the response variable is calculated. Second, the corresponding metric is evaluated using this predictive mean as point prediction. See here for a list of all supported metrics.

Tuning parameters aka hyperparameters for the tree boosting part

Below is a list of important parameters for the tree-boosting part. A comprehensive list of all tree-bosting related parameters can be found here.

• num_iterations 🔗︎, default = 100, type = int, aliases: num_iteration, n_iter, num_tree, num_trees, num_round, num_rounds, num_boost_round, n_estimators, constraints: num_iterations >= 0

  • number of boosting iterations

  • this is arguably the most important tuning parameter, in particular for regession settings

• learning_rate 🔗︎, default = 0.1, type = double, aliases: shrinkage_rate, eta, constraints: learning_rate > 0.0

  • shrinkage rate or damping parameter

  • smaller values lead to higher predictive accuracy but require more computational time since more boosting iterations are needed

• max_depth 🔗︎, default = -1, type = int

  • maximal depth of a tree

  • <= 0 means no limit

• num_leaves 🔗︎, default = 31, type = int, aliases: num_leaf, max_leaves max_leaf, constraints: 1 < num_leaves <= 131072

  • maximal number of leaves of a tree

• Note on ``max_depth`` and ``num_leaves`` parameters: The GPBoost library uses the LightGBM tree growing algorithm which grows trees using a leaf-wise strategy. I.e., trees are grown by first splitting leaf nodes that maximize the information gain until the maximal number of leaves num_leaves or the maximal depth of a tree max_depth is attained, even when this leads to unbalanced trees. This in contrast to a depth-wise growth strategy of other boosting implementations which builds “balanced” trees. For shallow trees (=small max_depth), there is likely no difference between these two tree growing strategies. If you only want to tune the maximal depth of a tree max_depth parameter and not the num_leaves parameter, it is recommended that you set the num_leaves parameter to a large value

• min_data_in_leaf 🔗︎, default = 20, type = int, aliases: min_data_per_leaf, min_data, min_child_samples, constraints: min_data_in_leaf >= 0

  • minimal number of samples in a leaf

• lambda_l2 🔗︎, default = 0.0, type = double, aliases: reg_lambda, lambda, constraints: lambda_l2 >= 0.0

  • L2 regularization

• lambda_l1 🔗︎, default = 0.0, type = double, aliases: reg_alpha, constraints: lambda_l1 >= 0.0

  • L1 regularization

• max_bin 🔗︎, default = 255, type = int, constraints: max_bin > 1

  • Maximal number of bins that feature values will be bucketed in

  • GPBoost uses histogram-based algorithms [1, 2, 3], which bucket continuous feature (covariate) values into discrete bins. A small number speeds up training and reduces memory usage but may reduce the accuracy of the model

• min_gain_to_split 🔗︎, default = 0.0, type = double, aliases: min_split_gain, constraints: min_gain_to_split >= 0.0

  • the minimal gain to perform a split

• line_search_step_length 🔗︎, default = false, type = bool

  • if true, a line search is done to find the optimal step length for every boosting update (see, e.g., Friedman 2001). This is then multiplied by the learning_rate

  • applies only to the GPBoost algorithm

• reuse_learning_rates_gp_model 🔗︎, default = true, type = bool

  • if true, the learning rates for the covariance and potential auxiliary parameters are kept at the values from the previous boosting iteration and not re-initialized when optimizing them

  • this option can only be used if optimizer_cov = gradient_descent or optimizer_cov = lbfgs (for the latter, the approximate Hessian is reused)

• train_gp_model_cov_pars 🔗︎, default = true, type = bool

  • if true, the covariance parameters of the Gaussian process / random effects model are trained (estimated) in every boosting iteration of the GPBoost algorithm, otherwise not

• use_gp_model_for_validation 🔗︎, default = true, type = bool

  • set this to true to also use the Gaussian process / random effects model (in addition to the tree model) for calculating predictions on the validation data when using the GPBoost algorithm

• leaves_newton_update 🔗︎, default = false, type = bool

  • if true, a Newton update step is done for the tree leaves after the gradient step

  • applies only to the GPBoost algorithm for Gaussian data and cannot be used for non-Gaussian data