---
documentclass: jss
classoption: article, shortnames, nojss
author:
  - name: Jiangyan Zhao 
    affiliation: |
      | School of Statistics
      | East China Normal University  
    address: |
      | 3663 North Zhongshan Road,
      | Shanghai 200062, China 
    email: \email{jyzhao@sfs.ecnu.edu.cn}
  - name: Kunhai Qing 
    address: |
      | 3663 North Zhongshan Road,
      | Shanghai 200062, China 
    email: \email{51254404017@stu.ecnu.edu.cn}
    affiliation: |
      | School of Statistics
      | East China Normal University \AND
    # To add another line, use \AND at the end of the previous one as above
  - name: Jin Xu 
    address: |
      | 3663 North Zhongshan Road,
      | Shanghai 200062, China 
    email: \email{jxu@stat.ecnu.edu.cn}
    affiliation: |
      | School of Statistics
      | \emph{and}
      | Key Laboratory of Advanced Theory and Application in Statistics and Data Science - MOE
      | East China Normal University 
title:
  formatted: "\\pkg{BKP}: An \\proglang{R} Package for Beta Kernel Process Modeling"
  # If you use tex in the formatted title, also supply version without
  plain:     "BKP: An R Package for Beta Kernel Process Modeling"
  # For running headers, if needed
  short:     "Beta Kernel Process Modeling in \\proglang{R}"
abstract: >
 We present \pkg{BKP}, a user-friendly and extensible \proglang{R} package that implements the Beta Kernel Process (BKP) -- a fully nonparametric and computationally efficient framework for modeling spatially varying binomial probabilities. The BKP model combines localized kernel-weighted likelihoods with conjugate beta priors, resulting in closed-form posterior inference without requiring latent variable augmentation or intensive MCMC sampling. The package supports binary and aggregated binomial responses, allows flexible choices of kernel functions and prior specification, and provides loss-based kernel hyperparameter tuning procedures. In addition, BKP extends naturally to the Dirichlet Kernel Process (DKP) for modeling spatially varying categorical or multinomial data. To our knowledge, this is the first publicly available \proglang{R} package for implementing BKP-based methods. We illustrate the use of \pkg{BKP} through several synthetic and real-world datasets, highlighting its interpretability, accuracy, and scalability. The package aims to facilitate practical application and future methodological development of kernel-based beta modeling in statistics and machine learning.
keywords:
  # at least one keyword must be supplied
  formatted: [Beta kernel process, Dirichlet kernel process, nonparametric Bayesian modeling, Beta-binomial, Dirichlet-multinomial, "\\proglang{R}"]
  plain:     [Beta kernel process, Dirichlet kernel process, nonparametric Bayesian modeling, Beta-binomial, Dirichlet-multinomial, R]
preamble: >
  \usepackage{amsmath}
  \usepackage{orcidlink,thumbpdf,lmodern} 
  \usepackage{framed}
  \usepackage{amsmath, amssymb, amsthm}
  \newtheorem{remark}{Remark}
  \usepackage{subcaption} 
  \usepackage{rotating, threeparttable, booktabs, caption, dcolumn, pdflscape}
  \usepackage{dsfont, pifont}
  \usepackage{soul}  
  \newcommand{\class}[1]{`\code{#1}'}
  \newcommand{\fct}[1]{\code{#1()}}
  \renewcommand{\vec}[1]{\mbox{\boldmath ${#1}$}}
  \def\argmin{\mbox{\rm argmin}}
  \def\argmax{\mbox{\rm argmax}}
output: 
  rticles::jss_article: 
    keep_tex: yes
bibliography: refs.bib
vignette: >
  %\VignetteIndexEntry{BKP: An R Package for Beta Kernel Process Modeling}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, setup, include=FALSE}
library(BKP)
library(tgp)
library(gplite)
library(kernlab)
library(mlbench)
library(pROC)
library(gridExtra)
library(mdhglm)
library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
options(prompt = "R> ", continue = "+  ", width = 70,
        useFancyQuotes = FALSE)
```

# Introduction {#sec:intro}
Estimating a continuous probability function from binary or binomial observations is a fundamental task in modern statistics and machine learning \citep{Hastie2009StatLearning, Rolland2018BKP, Murphy2022PML}. Such problems arise in various domains where the goal is to infer a latent probability surface from individual Bernoulli outcomes or aggregated binomial counts over a continuous input space. Representative applications include binary classification \citep{MacKenzie2014, Wen2025KLR}, probability calibration \citep{Sung2020binaryCalibration, Dimitriadis2023}, relative abundance modeling \citep{Martin2020}, and longitudinal analysis of patient-reported outcomes \citep{Najera2018HRQoL, Najera2019PRO}. 

Classical approaches such as logistic regression and generalized additive models provide computationally efficient tools for binomial regression, but their parametric form often limits their ability to capture complex nonlinear patterns \citep{Hastie2009StatLearning}. Nonparametric models like Gaussian Process (GP) classifiers provide greater modeling flexibility and principled uncertainty quantification \citep{Rasmussen2006GPML}, yet the non-Gaussian likelihood of binary responses necessitates approximate inference, leading to substantial computational cost and implementation complexity \citep{Nickisch2008BinaryGP}. 

The \pkg{BKP} package implements the \emph{Beta Kernel Process} (BKP), a scalable and interpretable nonparametric model for estimating binomial probability surfaces. Developed in \proglang{R} \citep{R} and available on the Comprehensive R Archive Network (CRAN) at \url{https://cran.r-project.org/package=BKP}, \pkg{BKP} combines localized kernel-weighted likelihoods with conjugate beta priors, yielding closed-form posterior inference without latent variables or numerical approximation. It supports both binary and aggregated binomial data, and delivers full posterior uncertainty quantification while maintaining linear scalability and computational efficiency \citep{Goetschalckx:2011, MacKenzie2014, Rolland2018BKP}.

Compared to GP-based methods, BKP eliminates the need for sampling or variational inference and provides interpretable, locally updated estimates. This makes it suitable for applications demanding transparency, real-time feedback, or adaptive decision-making. BKP has been extended to the multinomial setting via the \emph{Dirichlet Kernel Process} (DKP), enabling nonparametric modeling of categorical proportions across multiple groups. Despite its practical appeal, the theoretical properties of BKP remain underexplored. \citet{Mussi2024KMAB} identified the lack of formal guarantees for BKP-based decision processes as an open problem, emphasizing the need for further methodological development.

While GP models are supported by mature software libraries, such as \pkg{DiceKriging} \citep{Roustant2012DiceKriging}, \pkg{GPfit} \citep{MacDonald2015GPfit}, and \pkg{gplite} \citep{Piironen2022gplite} in \proglang{R}, and \pkg{GPyTorch} \citep{Gardner2018GPyTorch}, \pkg{GPflow} \citep{Matthews2017GPflow}, and \pkg{GPy} \citep{gpy2014} in \proglang{Python}. So far, there are about 100 GP-related packages on CRAN (as seen by running \code{packagefinder::findPackage("Gaussian process", display = "browser")}).  In contrast, software for kernel-based modeling of binomial and multinomial data remains scarce. The \pkg{BKP} package fills this gap by offering a natively, unified and extensible framework for nonparametric modeling of  binary/binomial, categorical/multinomial data. It provides multiple options for prior specification, kernel functions, and loss functions, enabling flexible model customization. These features make \pkg{BKP} applicable to a wide range of problems in biomedicine, ecology, social science, and industrial statistics. 

The remainder of the paper is organized as follows. Section \ref{sec:model} introduces the statistical foundation of the BKP model, including prior specification strategies, kernel functions, loss-based hyperparameter tuning, and the DKP extension. Section \ref{sec:package} describes the structure and main functionalities of the \pkg{BKP} package, including model fitting, prediction, and simulation for both BKP and DKP. Section \ref{sec:example} presents illustrative examples to demonstrate the use of BKP and DKP in binary classification and compositional modeling tasks. Section \ref{sec:Loaloa} demonstrates the BKP workflow using the Loa loa dataset. Section \ref{sec:summary} concludes with a discussion of current limitations and potential extensions.

# Statistical Foundation {#sec:model}
 
## Beta Kernel Process {#sec:review-bkp}

Let $\vec{x} = (x_1, x_2, \ldots, x_d) \in \mathcal{X} \subset \mathbb{R}^d$ denote a $d$-dimensional input. Suppose the success probability surface $\pi(\vec{x}) \in [0,1]$ is unknown. At each location $\vec{x}$, the observed data is modeled as
\[
	y(\vec{x}) \sim \mathrm{Binomial}(m(\vec{x}), \pi(\vec{x})),
\]
where $y(\vec{x})$ is the number of successes out of $m(\vec{x})$ independent trials. 
The full dataset comprises $n$ observations $\mathcal{D}_n = \{(\vec{x}_i, y_i, m_i)\}_{i=1}^n$, where we write $y_i= y(\vec{x}_i)$ and $m_i= m(\vec{x}_i)$ for brevity.

In line with the Bayesian paradigm, assign a Beta prior to the unknown probability function as
\[
	\pi(\vec{x}) \sim \mathrm{Beta}(\alpha_0(\vec{x}), \beta_0(\vec{x})),
\]
where $\alpha_0(\vec{x}) > 0$ and $\beta_0(\vec{x}) > 0$ are spatially varying shape parameters. Details on prior specification are discussed in Section \ref{sec:prior}.

Let $k: \mathcal{X} \times \mathcal{X} \to [0,1]$ denote a user-defined kernel function measuring the similarity between input locations. 
By the kernel-based Bayesian updating strategy \citep{Goetschalckx:2011, Rolland2018BKP}, the BKP model constructs a closed-form posterior distribution for $\pi(\vec{x})$ as
\begin{align}\label{eq:BKP_model}
	\pi(\vec{x}) \mid \mathcal{D}_n &\sim \mathrm{Beta}\left(\alpha_n(\vec{x}), \beta_n(\vec{x})\right), \notag \\
	\alpha_n(\vec{x}) &= \alpha_0(\vec{x}) + \sum_{i=1}^{n} k(\vec{x}, \vec{x}_i) y_i = \alpha_0(\vec{x}) + \vec{k}^\top(\vec{x}) \vec{y}, \\
	\beta_n(\vec{x})  &= \beta_0(\vec{x}) + \sum_{i=1}^{n} k(\vec{x}, \vec{x}_i) (m_i - y_i) = \beta_0(\vec{x}) + \vec{k}^\top(\vec{x}) (\vec{m} - \vec{y}), \notag
\end{align}
where $\vec{k}(\vec{x}) = [k(\vec{x}, \vec{x}_1), \ldots, k(\vec{x}, \vec{x}_n)]^\top$ is the vector of kernel weights and $\vec{y}=(y_1,\ldots,y_n)^\top$.
\begin{remark}
While the BKP model leverages the conjugacy of the Beta-Binomial pair, it differs from the traditional Bayesian paradigm in the sense that the posterior update is induced from a kernel-weighted local likelihood defined by
\[
	\widetilde{L}(\pi(\vec{x}); \mathcal{D}_n) 
	\propto \prod_{i=1}^{n} \left\{\pi(\vec{x})^{y_i}(1-\pi(\vec{x}))^{m_i - y_i} \right\}^{k(\vec{x}, \vec{x}_i)}= \pi(\vec{x})^{\vec{k}^\top(\vec{x}) \vec{y}} \{1 - \pi(\vec{x})\}^{\vec{k}^\top(\vec{x}) (\vec{m}-\vec{y})}.
\]
This approach mimics the local likelihood method of \citet{Fan1998LocalMLE}, where data are reweighted based on distance to the target point in the input space. 
And, the choice of kernel parameters is driven by empirical risk minimization rather than posterior inference. 
Thus, BKP is best interpreted as a nonparametric, Bayesian-inspired smoothing framework.
\end{remark}

\begin{remark}
One prominent advantage of the closed-form updating scheme is its light burden of computational complexity. Fitting the BKP model involves $\mathcal{O}(n^2)$ operations for computing the kernel matrix, in contrast to the $\mathcal{O}(n^3)$ operations typically required by Gaussian process regression. While, evaluating the posterior at a new location requires only $\mathcal{O}(n)$ kernel computations. 
\end{remark}


Based on the resulting posterior distribution (\ref{eq:BKP_model}), the posterior mean
\begin{equation}\label{eq:BKPmean}
	\widehat{\pi}_n(\vec{x}) = \mathbb{E}[\pi(\vec{x}) \mid \mathcal{D}_n] = \frac{\alpha_n(\vec{x})}{\alpha_n(\vec{x}) + \beta_n(\vec{x})}
\end{equation}
serves as a smooth estimator of the latent success probability. The corresponding posterior variance  
\begin{equation}\label{eq:BKPvar}
	\sigma^2_n(\vec{x})=\mathrm{Var}[\pi(\vec{x}) \mid \mathcal{D}_n] = \frac{\widehat{\pi}_n(\vec{x})\{1 - \widehat{\pi}_n(\vec{x})\}}{\alpha_n(\vec{x}) + \beta_n(\vec{x}) + 1}
\end{equation}
provides a local measure of epistemic uncertainty.
These posterior summaries can be used to visualize prediction quality across the input space, particularly highlighting regions with sparse data coverage. See Section \ref{sec:example} for illustrations.

For binary classification, the posterior mean can be thresholded to produce hard predictions through
\begin{equation}\label{eq:class_cri}
	\widehat{y}(\vec{x}) =
	\begin{cases}
		1 & \text{if } \widehat{\pi}_n(\vec{x}) > \pi_0, \\
		0 & \text{otherwise},
	\end{cases}
\end{equation}
where $\pi_0 \in (0,1)$ is a user-specified threshold, typically set to be 0.5.


## Prior Specification {#sec:prior}

We provide three strategies for prior specification as follows. 
\begin{itemize}
	\item \textbf{Non-informative prior:}  
	A default and widely used choice is the uniform prior, which sets $\alpha_0(\vec{x}) = \beta_0(\vec{x}) \equiv 1$ for all $\vec{x}$. %This prior is minimally informative and corresponds to a flat density over the interval $(0,1)$, reflecting complete prior ignorance. 
    It is appropriate when no prior knowledge is available.
	
	\item \textbf{Informative prior with fixed mean:}  
	When prior information about the overall success probability $p_0 \in (0,1)$ is available, an informative prior can be constructed by setting
	\[
		\alpha_0(\vec{x}) = r_0 p_0, \qquad \beta_0(\vec{x}) = r_0 (1 - p_0),
	\]
	where $r_0 > 0$ is a scalar precision parameter controlling the strength of the prior. Larger value of $r_0$ represents greater prior certainty. It contains the previous non-informative prior as a special case with $r_0 = 2$ and $p_0 = 0.5$. In practice, we recommend specifying $p_0$ according to the empirical overall success proportion $y_i/m_i$.
	
	\item \textbf{Data-adaptive informative prior:}  
	To accommodate spatial variation in the data, \pkg{BKP} supports a locally adaptive prior in which both the prior mean and prior strength vary across the input space. Specifically, at each location $\vec{x}$, the prior mean $p(\vec{x})$ and prior precision $r(\vec{x})$ are estimated using kernel-weighted averages through
	\[
		p(\vec{x}) = \sum_{i=1}^n w_i(\vec{x}) \cdot \frac{y_i}{m_i}, \qquad r(\vec{x}) = r_0 \sum_{i=1}^n k(\vec{x}, \vec{x}_i),
	\]
	where $r_0 > 0$ is a global precision parameter and $w_i(\vec{x})$ are normalized kernel weights defined by $w_i(\vec{x}) = {k(\vec{x}, \vec{x}_i)}/{\sum_{j=1}^n k(\vec{x}, \vec{x}_j)}$.
	Then, the prior parameters at $\vec{x}$ are given by
	\[
	\alpha_0(\vec{x}) = r(\vec{x}) p(\vec{x}), \qquad \beta_0(\vec{x}) = r(\vec{x}) \{1 - p(\vec{x})\}.
	\]
	
	This data-adaptive formulation allows the prior to dynamically respond to the local sampling density. In well-sampled regions, it becomes more concentrated, while in sparse regions, it remains diffuse. Such adaptivity improves calibration in under-sampled areas and reduces overfitting where data are dense, thereby enhancing inference robustness under spatial heterogeneity. This strategy is consistent with the principle of local likelihood modeling \citep{Fan1996LocalModel}.
\end{itemize}

We offer some practical guideline to choose the global precision parameter. 
In classification when each input location receives a single categorical observation, we recommend choosing a relatively small value, such as $0.01 \leq r_0 \leq 0.1$, to prevent the prior domination and keep the posterior inference primarily data-driven. We illustrate this point in Example~1 of Section \ref{sec:example_bkp}. 
For model fitting with binomial responses at each location, a reasonable choice is to set $r_0$ to the mean number of trials per location, reflecting the typical information content of the data. 
In the absence of reliable prior information, a simple default is $r_0 = 0.1$, with subsequent adjustment guided by cross-validation or predictive performance over a range of candidate values (e.g., $r_0 \in \{0.001, 0.01, 0.1, 1, 5, 10\}$) to assess the robustness of predictions and decision boundaries.

## Model Selection via Kernel Hyperparameter Tuning {#sec:model-selection}

### Kernel Functions 

Let $h(\vec{x}, \vec{x}'; \vec{\theta})$ denote the scaled Euclidean distance:
\[
h(\vec{x}, \vec{x}'; \vec{\theta}) = \sqrt{\sum_{j=1}^{d} \left(\frac{x_{j} - x_{j}'}{\theta_{j}}\right)^{2}},
\]
where $\vec{\theta} = (\theta_1, \theta_2, \ldots, \theta_d)$ are positive kernel hyperparameters governing the relative importance of each input component. Based on this metric, define the kernel function as $k(\vec{x}, \vec{x}') = k(h)$, where the functional form of $k(\cdot)$ determines the kernel type \citep{Rasmussen2006GPML}.

The \pkg{BKP} package implements several widely used kernel functions, namely Gaussian (squared-exponential) and Matérn families, summarized in Table \ref{tab:kernel-functions}.

\begin{table}[!t]
	\renewcommand{\arraystretch}{1.5}
	\centering
	\caption{Kernel functions implemented in \pkg{BKP}.}
	\label{tab:kernel-functions}
	\begin{tabular}{ll}
		\toprule
		\textbf{Kernel Type} & \textbf{Function $k(h)$}\\
		\midrule
		Gaussian & $k(h) = \exp(-h^2)$ \\
		Mat\'{e}rn $\nu = 5/2$ & $k(h) = \left(1 + \sqrt{5}h + \frac{5}{3}h^2 \right)\exp(-\sqrt{5}h)$ \\
		Mat\'{e}rn $\nu = 3/2$ & $k(h) = \left(1 + \sqrt{3}h\right)\exp(-\sqrt{3}h)$ \\
		\bottomrule
	\end{tabular}
\end{table}

### Loss Functions  

Hyperparameter tuning is conducted by minimizing a user-specified loss function evaluated using the leave-one-out cross-validation (LOOCV) method \citep{Montesano2012}. For each data point $\vec{x}_i$, the model is refitted up on the remaining data $\mathcal{D}_n^{-i}$, and the posterior mean $\widehat{\pi}_n^{-i}(\vec{x}_i;\vec{\theta})$ is used as the prediction. The LOOCV procedure is known to mitigate overfitting and produce better generalization performance than marginal likelihood-based approaches \citep{Rasmussen2006GPML, Vehtari2017LOOCV}.

Currently, \pkg{BKP} supports two loss functions: the Brier score (i.e., squared error loss) and the log-loss (i.e., negative log-likelihood or cross-entropy).

\paragraph{Brier Score}
Let $\widetilde{\pi}_i = y_i / m_i$ denote the empirical success proportion at input location $\vec{x}_i$. 
The LOOCV Brier score is defined as
\begin{equation}\label{eq:Brier}
	\mathrm{BS}(\vec{\theta}; \mathcal{D}_n) = \frac{1}{n} \sum_{i=1}^{n} \left\{ \widehat{\pi}_n^{-i}(\vec{x}_i;\vec{\theta}) - \widetilde{\pi}_i \right\}^2,
\end{equation}
which penalizes the squared deviation between the predicted and observed success proportions.

\paragraph{Log-Loss}
The log-loss is widely used for probabilistic classification. Based on the LOOCV predictions, it is defined as
\begin{equation}\label{eq:log-loss}
	\mathrm{LL}(\vec{\theta}; \mathcal{D}_n) = -\frac{1}{n} \sum_{i=1}^{n} \left[ \widetilde{\pi}_i \log \widehat{\pi}_n^{-i}(\vec{x}_i;\vec{\theta}) + (1 - \widetilde{\pi}_i) \log \left\{1 - \widehat{\pi}_n^{-i}(\vec{x}_i;\vec{\theta})\right\} \right].
\end{equation}
It corresponds to the negative log-likelihood of the binomial model evaluated at the LOOCV predictive probabilities.

\begin{remark} %{Comparison of the Two Criteria.}
Although both criteria rely on LOOCV estimates, they emphasize different aspects of predictive performance \citep{Gneiting2007ProperScoring, Flores2025ClassifierEvaluation}. 
The Brier score directly penalizes squared deviations from empirical proportions, making it more robust to poorly calibrated probabilities and better suited for smooth probability estimation under binomial data. 
In contrast, the log-loss penalizes overconfident mispredictions more heavily and is highly sensitive to the accuracy of the predicted probabilities. In practice, the choice between the two depends on whether calibration quality or robustness is prioritized in the evaluation \citep{DRatings2023}. 
\end{remark}

### Hyperparameter Optimization {#sec:HPO}

To improve optimization stability and reduce irregularities near the boundary of the parameter space, we adopt the reparameterization strategy proposed by \citet{MacDonald2015GPfit}. Specifically, we transform the kernel scale parameters via $\gamma_j = \log_{10}(\theta_j)$ for $j = 1, \ldots, d$, and denote $\vec\gamma = (\gamma_1, \ldots, \gamma_d)$. For example, the Gaussian kernel becomes
\begin{equation}\label{eq:reparameterization}
	k(\vec{x}, \vec{x}'; \vec\gamma) = \exp\left\{ -\sum_{j=1}^{d}\left(\frac{x_j - x_j'}{10^{\gamma_j}} \right)^2 \right\}.
\end{equation}

Although gradient-based optimizers such as L-BFGS-B \citep{Byrd:1995} are computationally efficient, their performance can be sensitive to initialization in non-convex settings. 
To improve robustness, we adopt a multi-start strategy based on a Latin Hypercube design (LHD) of $n_0 = 10d$ initial runs within a carefully constructed region \citep{Loeppky:2009}.

By \citet{MacDonald2015GPfit}, the Gaussian kernel value $k(\vec{x}, \vec{x}'; \vec\gamma)$ typically lies within the interval $[0.0067, 0.9999]$, approximately $[\exp(-5), \exp(-10^{-4})]$, to ensure stable numerical behavior. Assuming a space-filling design of size $n = 10d$ over $[0,1]^d$ and isotropic smoothness (i.e., $\gamma_j = \gamma_0$ for all $j$),the minimum pairwise distance along each coordinate is roughly $|x_{ik} - x_{jk}| \approx 0.1$. 
This leads to the following search region for $\vec\gamma$:
\begin{equation}\label{eq:gamma-bound}
	\Omega_0 = \left[\frac{\log_{10} d - \log_{10} 500}{2},\; \frac{\log_{10} d + 2}{2} \right]^d.
\end{equation}
For the mat\'{e}rn kernels, we adopt the same search region.

The initial values ${\vec{\gamma}^{(1)}, \ldots, \vec{\gamma}^{(n_0)}} \subset \Omega_0$ are drawn using LHD. For each initial value, solve
\[
	\widehat{\vec{\gamma}}^{(i)} = \underset{\vec{\gamma} \in \Omega}{\argmin}\; L(\vec{\gamma};\mathcal{D}_n), \quad i=1,\ldots, n_0,
\]
where $L(\cdot)$ is the chosen loss function from \eqref{eq:Brier} or \eqref{eq:log-loss}, and $\Omega = [-3, 3]^d$  is a broader search region than $\Omega_0$. The final estimate is chosen as the best local minimum through
\[
	\widehat{\vec{\gamma}} = \underset{1 \leq i \leq n_0}{\argmin}\;  L(\widehat{\vec{\gamma}}^{(i)};\mathcal{D}_n).
\]

The procedure is summarized below:
\begin{enumerate}
	\item Generate $10d$ initial values $\vec{\gamma}$ in $\Omega_0$ using a space-filling LHD.
	\item Run the L-BFGS-B algorithm from each initial value.
	\item Select the solution that produces the lowest loss function.
\end{enumerate}

\begin{remark}
The implementation of LOOCV necessitates refitting the model $n$ times, which can be computationally intensive. However, the closed-form posterior updates in the BKP model mitigate this cost by reducing the computational complexity to $\mathcal{O}(n^2)$. This represents a substantial improvement over the typical $\mathcal{O}(n^3)$ complexity of Gaussian process models with binomial likelihoods, as demonstrated in Example 3.
\end{remark}


## Extension to Dirichlet Kernel Process {#sec:dkp}

The BKP model can be naturally extended to handle multi-class responses via the Dirichlet Kernel Process (DKP) on replacing the binomial likelihood with a multinomial model and the Beta prior with a Dirichlet prior \citep{MacKenzie2014}. 

Let the response at input $\vec{x} \in \mathcal{X} \subset \mathbb{R}^d$ be $\vec{y}(\vec{x}) = \left(y_1(\vec{x}), \ldots, y_q(\vec{x})\right)$, where $y_s(\vec{x})$ denotes the count of class $s$ out of $m(\vec{x}) = \sum_{s=1}^q y_s(\vec{x})$ total trials. Assume
\[
	\vec{y}(\vec{x}) \sim \mathrm{Multinomial}(m(\vec{x}), \vec{\pi}(\vec{x})),
\]
with class probabilities $\vec{\pi}(\vec{x}) = (\pi_1(\vec{x}), \ldots, \pi_q(\vec{x}))$ satisfying $\sum_{s=1}^q \pi_s(\vec{x}) = 1$.

A Dirichlet prior is imposed on $\vec{\pi}(\vec{x})$ as
\[
	\vec{\pi}(\vec{x}) \sim \mathrm{Dirichlet}(\vec{\alpha}_0(\vec{x})),
\]
where $\vec{\alpha}_0(\vec{x}) = (\alpha_{0,1}(\vec{x}), \ldots, \alpha_{0,q}(\vec{x}))$ are prior concentration parameters. Given training data $\mathcal{D}_n = \{(\vec{x}_i, \vec{y}_i)\}_{i=1}^n$, define the response matrix as $\vec{Y} = [\vec{y}_1, \ldots, \vec{y}_n]^\top \in \mathbb{R}^{n \times q}$. Then the kernel-smoothed conjugate posterior becomes
\begin{equation}\label{eq:DKP_model}
	\vec{\pi}(\vec{x}) \mid \mathcal{D}_n \sim \mathrm{Dirichlet}\left(\vec{\alpha}_n(\vec{x})\right), \quad \text{with} \quad \vec{\alpha}_n(\vec{x}) = \vec{\alpha}_0(\vec{x}) + \vec{k}^\top(\vec{x}) \vec{Y}.
\end{equation}

The posterior means and variances of the class probabilities are
\[
	\widehat{\pi}_{n,s}(\vec{x}) = \frac{\alpha_{n,s}(\vec{x})}{\sum_{s'=1}^q \alpha_{n,s'}(\vec{x})}, \quad \sigma^2_{n,s}(\vec{x}) = \frac{\widehat{\pi}_{n,s}(\vec{x})\{1 - \widehat{\pi}_{n,s}(\vec{x})\}}{\sum_{s'=1}^q \alpha_{n,s'}(\vec{x}) + 1}, \quad s = 1, \ldots, q.
\]
For classification tasks, labels are assigned by the maximum a posteriori (MAP) decision rule through
\begin{equation}\label{eq:class_cri_DKP}
	\widehat{y}(\vec{x}) = \underset{s \in \{1,\ldots,q\}}{\argmax}\; \widehat{\pi}_{n,s}(\vec{x}).
\end{equation}

Prior choices (e.g., non-informative, fixed informative, or locally adaptive) follow similarly from the BKP framework in Section \ref{sec:prior} by treating component-wise specification of prior class proportions.

Kernel hyperparameters are tuned by minimizing LOOCV-based loss functions, as described in Section \ref{sec:model-selection}. For multi-class responses, we consider the Brier score and the log-loss, defined respectively as
\begin{align}
	\mathrm{BS}_{\mathrm{multi}}(\vec{\theta}) &= \frac{1}{n} \sum_{i=1}^n \sum_{s=1}^q \left( \widehat{\pi}_{n,s}^{-i}(\vec{x}_i;\vec{\theta}) - \widetilde{\pi}_{i,s}\right)^2, \\
	\mathrm{LL}_{\mathrm{multi}}(\vec{\theta}) &= -\frac{1}{n} \sum_{i=1}^n \sum_{s=1}^q \widetilde{\pi}_{i,s} \log \widehat{\pi}_{n,s}^{-i}(\vec{x}_i;\vec{\theta}),
\end{align}
where $\widetilde{\pi}_{i,s}=y_{i,s}/m_i$ denotes the empirical class proportion.
These multi-class criteria reduce to their binary counterparts when $q=2$. %For the log-loss, the expressions coincide exactly, while for the Brier score one obtains $\mathrm{BS}{\mathrm{multi}}(\vec{\theta})\big|{q=2} = 2 \cdot \mathrm{BS}(\vec{\theta})$, due to the summation over both classes.

# BKP Package {#sec:package}

The \pkg{BKP} package provides the implementation of the Beta Kernel Process (BKP) and its extension to the Dirichlet Kernel Process (DKP). The package is designed with usability and extensibility in mind, offering a modular structure with separate functions for model fitting, prediction, simulation, kernel construction, prior specification, and hyperparameter tuning. Users can specify the kernel type and prior form to adapt to various application scenarios.
The core functionality of the package is summarized in Table \ref{tab:main_functions}.

The package can be obtained from CRAN using \code{install.packages("BKP")} or the development version from GitHub via \code{pak::pak("Jiangyan-Zhao/BKP")}.

\begin{table}[!t]
	\centering
	\caption{Main functions and S3 methods provided in the \pkg{BKP} package.}
	\label{tab:main_functions}
	\begin{tabular}{ll}
		\toprule
		 & Description \\
		\midrule
		\multicolumn{2}{l}{\textbf{Functions}} \\
		\fct{fit\_BKP} & Fit a Beta Kernel Process (BKP) model for binary/binomial data. \\
		\fct{fit\_DKP} & Fit a Dirichlet Kernel Process (DKP) model for categorical or multinomial data. \\
		\midrule
		\multicolumn{2}{l}{\textbf{S3 Methods}} \\
		\fct{predict} & Posterior prediction for BKP and DKP models. \\
		\fct{simulate} & Posterior simulation for BKP and DKP models. \\
		\fct{summary} & Summarize fitted BKP and DKP models. \\
		\fct{plot} & Visualization of fitted BKP and DKP models. \\
		\fct{print} & Print fitted models, predictions, simulations, and summaries. \\
		\bottomrule
	\end{tabular}
\end{table}


## BKP Model 

We first provide a detailed description of the BKP model. The DKP model, presented in Section \ref{sec:DKP_pkg}, highlights only the aspects that differ from the BKP implementation.

The function \fct{fit\_BKP} fits the BKP model \eqref{eq:BKP_model} for binomial or binary response data. Its main arguments are
\begin{Code}
fit_BKP(X, y, m, Xbounds = NULL, 
	prior = "noninformative", r0 = 2, p0 = mean(y / m),
	loss = "brier",  kernel = "gaussian", n_multi_start = NULL, theta = NULL)
\end{Code}

The first group specifies the core data inputs: \code{X} is an $n \times d$ input matrix, \code{y} is the vector of observed successes, and \code{m} is the corresponding vector of total trials. 
\code{Xbounds} rescales \code{X} to the unit hypercube $[0,1]^d$ for kernel hyperparameter optimization, and if \code{NULL} (default), \code{X} is assumed standardized.

The second group defines the prior for the latent success probability surface. \code{prior} selects the type (\code{"noninformative"}, \code{"fixed"}, or \code{"adaptive"}; see Section \ref{sec:prior}), and the prior is constructed using the \fct{get\_prior} function. \code{r0} and \code{p0} specify prior precision and mean when applicable.

The third group specifies the loss function and optimization settings. \code{loss} chooses the criterion for kernel hyperparameter estimation (\code{"brier"} (default) or \code{"log_loss"}), which is implemented via the \fct{loss\_fun} function, and \code{kernel} selects the kernel function (\code{"gaussian"} (default), \code{"matern32"}, or \code{"matern52"}). \code{n_multi_start} sets the number of random initial points for multi-start optimization (default $10d$). 
Alternatively, \code{theta} can fix the kernel length scale directly, as a scalar (applied to all dimensions) or a vector of length $d$, skipping optimization.
 Optimization is performed using \fct{multistart} from the \pkg{optimx} package \citep{Nash2011optimx, Nash2014optim}.
Note that both \fct{get\_prior} and \fct{loss\_fun} are primarily invoked internally by \fct{fit\_BKP}, but exporting them allows advanced users to inspect the construction of priors or customize the evaluation of the loss criterion if desired.

\begin{remark}
The current version is implemented entirely in \proglang{R}. Computational efficiency can be further improved by implementing additional components in \proglang{C++} via \pkg{Rcpp} \citep{Eddelbuettel2011Rcpp}, and by enabling options that leverage parallel computing architectures, for example, to accelerate multi-start optimization routines used in hyperparameter tuning.
\end{remark}

\fct{fit\_BKP} returns an object of class \class{BKP}, containing the estimated model hyperparameters $\vec{\theta}$ (\code{theta_opt}); the minimum achieved loss value (\code{loss_min}); the prior parameters $\alpha_{0}(\vec{x})$ and $\beta_{0}(\vec{x})$ (\code{alpha_0} and \code{beta_0}); the fitted posterior parameters $\alpha_{n}(\vec{x})$ and $\beta_{n}(\vec{x})$ (\code{alpha_n} and \code{beta_n}); and other input arguments.

The resulting \class{BKP} object can be directly used as the \code{object} argument in the functions \fct{predict}, \fct{simulate}, and \fct{summary}, or as the \code{x} argument in the functions \fct{plot} and \fct{print}.

For a \class{BKP} object \code{BKPmodel}, the function call
\begin{Code}
predict(BKPmodel, Xnew = NULL, CI_level = 0.95, threshold = 0.5, ...)
\end{Code}
returns, for each $\vec{x}$ in \code{Xnew} (or in the training data if \code{Xnew = NULL}), the predictive mean $\widehat{\pi}_n(\vec{x})$, the predictive variance $\sigma^2_n(\vec{x})$, the lower and upper bounds of the \code{CI_level} credible interval, and the predictive parameters $\alpha_{n}(\vec{x})$ and $\beta_{n}(\vec{x})$. 
If \code{m = 1}, it also returns the predictive binary label (0 or 1), determined by comparing the posterior mean to the specified \code{threshold} according to the classification rule in \eqref{eq:class_cri}.

The function call
\begin{Code}
simulate(BKPmodel, nsim = 1, seed = NULL, Xnew = NULL, threshold = NULL, ...)
\end{Code}
generates random samples from the posterior predictive distribution of a fitted BKP model at the input locations in \code{Xnew} (or at the training data if \code{Xnew = NULL}).
Specifically, it draws \code{nsim} samples from the posterior Beta distribution at each row of \code{Xnew}, representing the distribution of success probabilities implied by the fitted model. 
If \code{threshold} is specified and \code{m = 1}, binary class labels (0 or 1) are generated by thresholding the simulated probabilities. 
An optional \code{seed} argument ensures reproducibility. 
Such posterior samples are useful for downstream decision-making tasks, for example Bayesian optimization via Thompson sampling \citep{Garnett2023BO}.

The function call
\begin{Code}
summary(BKPmodel, ...)
\end{Code}
produces a structured summary of a fitted \class{BKP} object. 
It reports the number of training observations, input dimensionality, kernel type, prior specification, and the optimized kernel hyperparameters, together with the minimum achieved loss value. 
In addition, it provides posterior summaries at the training points, including the posterior mean and variance of the success probabilities. 
This function offers a convenient way to inspect the fitted model configuration and key inferential results.

To facilitate user interaction, the core \class{BKP} object and the outputs of \fct{predict}, \fct{simulate}, and \fct{summary} are all implemented as S3 class objects. Specifically, the \class{BKP} object itself is equipped with a customized \fct{print} method for concise display, while the outputs of \fct{predict}, \fct{simulate}, and \fct{summary} are returned as dedicated classes (i.e., \class{predict\_BKP}, \class{simulate\_BKP}, \class{summary\_BKP}). Assigning these classes not only provides customized \fct{print} methods for human-readable summaries, but also supports future extensibility and backward-compatible addition of new S3 methods, such as \fct{plot} or specialized extraction functions.

The \fct{plot} method provides graphical visualization of a fitted \class{BKP} object, adapting to the input dimensionality. For one-dimensional inputs ($d=1$), it displays the posterior mean of the predicted success probability as a line plot, with a shaded region representing the 95\% credible interval and the observed proportions ($\widetilde{\pi}_i = y_i / m_i$) overlaid as points; for classification tasks, an optional curve of the maximum posterior probability can be added to illustrate decision confidence. For two-dimensional inputs ($d=2$), the method generates contour plots over a 2D prediction grid, by default producing four summary surfaces: the predictive mean, the predictive variance, and the 2.5\% and 97.5\% quantiles of the posterior distribution, providing a comprehensive view of predicted probabilities and associated uncertainty; if the argument \code{only_mean = TRUE} is specified, only the contour plot of the predictive mean surface is displayed. For high-dimensional inputs ($d>2$), users must specify which dimensions to plot via the \code{dims} argument.

In addition to the core methods in Table \ref{tab:main_functions}, the \pkg{BKP} package provides auxiliary methods for inspecting fitted models, including \fct{fitted} to extract fitted values (i.e., predictive means), \fct{parameter} to retrieve model parameters ($\alpha_{n}(\vec{x})$ and $\beta_{n}(\vec{x})$), and \fct{quantile} to obtain posterior quantiles, with the quantile levels specified by the argument \code{probs}, from a fitted \class{BKP} object. These functions complement the main workflow and facilitate exploration of model outputs and posterior summaries.

## DKP Model {#sec:DKP_pkg}

The function \fct{fit\_DKP} implements the DKP model~\eqref{eq:DKP_model}, extending the BKP framework to multinomial response data. The response \code{Y} should be an $n \times q$ matrix of observed counts, where each row represents an observation and each column a category. The \code{m} argument is not required, as the total count is implicitly defined by the row sums of \code{Y}. When using a fixed prior (\code{prior = "fixed"}), the user must provide a global prior mean vector \code{p0} of length $q$, which defaults to the empirical class proportions \code{colMeans(Y / rowSums(Y))}. 

The function returns an object of class \class{DKP}, which contains elements analogous to those in a \class{BKP} object, including optimized kernel hyperparameters, the selected kernel function, the loss specification, and prior information. Distinctively, the DKP posterior parameters are vectors $\vec{\alpha}_n(\vec{x})$ associated with each input $\vec{x}$, representing the smoothed Dirichlet parameters across categories. Other elements, such as the original and normalized input matrices, normalization bounds, and observed counts, are also included to support downstream prediction, simulation, and visualization.

Prediction, simulation, and visualization methods extend directly from the class \class{BKP}. The \fct{predict} method returns per-class predictive means, variances, and credible intervals; if all row sums of \code{Y} are 1, it additionally provides predicted class labels following the MAP classification rule in \eqref{eq:class_cri_DKP}. The \fct{simulate} method draws samples from the posterior predictive Dirichlet distribution and can optionally generate sampled class labels via Thompson sampling. For one-dimensional inputs ($d = 1$), \fct{plot} displays one curve per class with credible intervals, while for two-dimensional inputs ($d = 2$), it produces a panel of contour plots illustrating the predictive mean or other posterior summaries for each class.


# Examples using BKP {#sec:example}

```{r, Examples using BKP, include=FALSE}
# ============================================================== #
# ========================= BKP Examples ======================= #
# ============================================================== #
#-------------------------- 1D Example ---------------------------
#-------------------------- Example 1 ----------------------------
set.seed(123)
# Define true success probability function
true_pi_fun <- function(x) {
  1/(1+exp(-3*x))
}
# Data points
n <- 7
Xbounds <- matrix(c(-2,2), nrow = 1)
X <- lhs(n = n, rect = Xbounds)
true_pi <- true_pi_fun(X)
m <- sample(100, n, replace = TRUE)
y <- rbinom(n, size = m, prob = true_pi)

# Fit BKP model
BKP_model_1D_1 <- fit_BKP(X, y, m, Xbounds = Xbounds)

# Print BKP model
print(BKP_model_1D_1)

# New data points
Xnew = matrix(seq(-2, 2, length = 100), ncol = 1)
true_pi <- true_pi_fun(Xnew)

# Plot results
pdf("ex1.pdf", width = 6, height = 6)
plot(BKP_model_1D_1)
lines(Xnew,true_pi, col = "black", lwd = 2)
legend(x = -2.24,y = 0.89,
       legend = "True Probability",
       col = "black",
       lwd = 2,
       bty = "n",
       inset = 0.02)
dev.off()

# Simulate
sim <- simulate(BKP_model_1D_1, nsim = 3, Xnew = Xnew)
pdf("ex1_sim.pdf", width = 6, height = 6)
matplot(Xnew, sim$samples, type = "l", lty = 1, lwd = 1.5,
        col = rainbow(ncol(sim$samples)),
        xlab = "x", ylab = "Probability",
        main = "Simulated Probability Curves")
legend("topleft", legend = paste("sample", 1:ncol(sim$samples)),
       col = rainbow(ncol(sim$samples)), lty = 1, lwd = 2, bty = "n")
dev.off()

### Classification
set.seed(123)
# Data points
n <- 20
Xbounds <- matrix(c(-2,2), nrow = 1)
X <- lhs(n = n, rect = Xbounds)
true_pi <- true_pi_fun(X)
m <- rep(1,n)
y <- as.numeric(true_pi > 0.5)

# Fit BKP model with r0 = 0.01
BKP_model_1D_1_class_1 <- fit_BKP(
  X, y, m, Xbounds = Xbounds,
  prior = "fixed", r0 = 0.01, loss = "log_loss")

# New data points
Xnew = matrix(seq(-2, 2, length = 100), ncol = 1)
true_pi <- true_pi_fun(Xnew)

# Plot results
pdf("ex1class001.pdf", width = 6, height = 6)
plot(BKP_model_1D_1_class_1)
lines(Xnew,true_pi, col = "black", lwd = 2)
dev.off()

# Fit BKP model with r0 = 2
BKP_model_1D_1_class_2 <- fit_BKP(
  X, y, m, Xbounds = Xbounds,
  prior = "fixed", r0 = 2, loss = "log_loss")

# Plot results
pdf("ex1class2.pdf", width = 6, height = 6)
plot(BKP_model_1D_1_class_2)
lines(Xnew,true_pi, col = "black", lwd = 2)
dev.off()

#-------------------------- Example 2 ----------------------------
set.seed(123)
# Define true success probability function
true_pi_fun <- function(x) {
  (1 + exp(-x^2) * cos(10 * (1 - exp(-x)) / (1 + exp(-x)))) / 2
}
# Data points
n <- 30
Xbounds <- matrix(c(-2,2), nrow = 1)
X <- lhs(n = n, rect = Xbounds)
true_pi <- true_pi_fun(X)
m <- sample(100, n, replace = TRUE)
y <- rbinom(n, size = m, prob = true_pi)

# Fit BKP model
BKP_model_1D_2 <- fit_BKP(X, y, m, Xbounds = Xbounds)

# New data points
Xnew = matrix(seq(-2, 2, length = 100), ncol=1)
true_pi <- true_pi_fun(Xnew)

# Plot results
pdf("ex2.pdf", width = 6, height = 6)
plot(BKP_model_1D_2)
lines(Xnew,true_pi, col = "black", lwd = 2)
legend(x = -2.24,y = 0.89,
       legend = "True Probability",
       col = "black",
       lwd = 2,
       bty = "n",
       inset = 0.02)
dev.off()

# Simulate
sim <- simulate(BKP_model_1D_2, nsim = 3, Xnew = Xnew)
pdf("ex2_sim.pdf", width = 6, height = 6)
matplot(Xnew, sim$samples, type = "l", lty = 1, lwd = 1.5,
        col = rainbow(ncol(sim$samples)),
        xlab = "x", ylab = "Probability",
        main = "Simulated Probability Curves")
legend("topleft", legend = paste("sample", 1:ncol(sim$samples)),
       col = rainbow(ncol(sim$samples)), lty = 1, lwd = 2, bty = "n")
dev.off()

#-------------------------- 2D Example ---------------------------
#-------------------------- Example 3 ----------------------------
set.seed(123)
# Define 2D latent function and probability transformation
true_pi_fun <- function(X) {
  if(is.null(nrow(X))) X <- matrix(X, nrow=1)
  m <- 8.6928
  s <- 2.4269
  x1 <- 4*X[,1]- 2
  x2 <- 4*X[,2]- 2
  a <- 1 + (x1 + x2 + 1)^2 *
    (19- 14*x1 + 3*x1^2- 14*x2 + 6*x1*x2 + 3*x2^2)
  b <- 30 + (2*x1- 3*x2)^2 *
    (18- 32*x1 + 12*x1^2 + 48*x2- 36*x1*x2 + 27*x2^2)
  f <- log(a*b)
  f <- (f- m)/s
  return(pnorm(f))  # Transform to probability
}
# Data points
n <- 100
Xbounds <- matrix(c(0, 0, 1, 1), nrow = 2)
X <- lhs(n = n, rect = Xbounds)
true_pi <- true_pi_fun(X)
m <- sample(100, n, replace = TRUE)
y <- rbinom(n, size = m, prob = true_pi)

# Fit BKP model
BKP_model_2D <- fit_BKP(X, y, m, Xbounds=Xbounds)

# Print BKP model
print(BKP_model_2D)

# Plot results
pdf("ex3.pdf", width = 9, height = 8)
plot(BKP_model_2D)
dev.off()

# Plot True distribution
pdf("ex3_true.pdf", width = 4.5, height = 4)
Xnew1 <- seq(Xbounds[1,1], Xbounds[1,2], length.out = 200)
Xnew2 <- seq(Xbounds[2,1], Xbounds[2,2], length.out = 200)
Xnew <- expand.grid(Xnew1 = Xnew1, Xnew2 = Xnew2)
true_pi <- true_pi_fun(Xnew)
df <- data.frame(x1 = Xnew$Xnew1, x2 = Xnew$Xnew2,
                 True = true_pi)
print(BKP:::my_2D_plot_fun("True", title = "True Probability", data = df))
dev.off()

# Predict via BKP model
Xnew <- lhs(n = 10, rect = Xbounds)
predict(BKP_model_2D, Xnew)


# ============================================================
# Execution mode:
# run_elapsed_time = FALSE means using the pre-computed average
# elapsed time matrix, without re-running the time measurement
# experiments.
#
# Setting run_elapsed_time = TRUE will run the full benchmarking
# procedure for 5 methods × 5 n values × 20 repetitions.
# This can take many hours (potentially tens of hours) depending
# on your machine's performance.
#
# It is recommended to keep this FALSE unless you need to
# fully reproduce the timing experiments.
# ============================================================

run_elapsed_time = FALSE
if(run_elapsed_time){
  n_vals <- c(200, 500, 1000, 2000, 5000)
  elapsed_time_given <- matrix(NA, nrow = 20, ncol = 5)
  elapsed_time_single_start <- matrix(NA, nrow = 20, ncol = 5)
  elapsed_time_multi_start <- matrix(NA, nrow = 20, ncol = 5)
  elapsed_time_given_gp <- matrix(NA, nrow = 20, ncol = 5)
  elapsed_time_optim_gp <- matrix(NA, nrow = 20, ncol = 5)

  for (i in 1:20) {
    for (j in 1:5) {
      set.seed((i - 1) * 5 + j)
      cat("i =", i, "j =", j, ":")
      n <- n_vals[j]
      Xbounds <- matrix(c(0, 0, 1, 1), nrow = 2)
      X <- lhs(n = n, rect = Xbounds)
      true_pi <- true_pi_fun(X)
      m <- sample(100, n, replace = TRUE)
      y <- rbinom(n, size = m, prob = true_pi)

      elapsed_time_given[i, j] <- system.time({
        fit <- fit_BKP(X, y, m, Xbounds = Xbounds, theta = 1)
      })["elapsed"]
      cat(" BKP_given =", elapsed_time_given[i, j])

      elapsed_time_single_start[i, j] <- system.time({
        fit <- fit_BKP(X, y, m, Xbounds = Xbounds, n_multi_start = 1)
      })["elapsed"]
      cat(", BKP_single_start =", elapsed_time_single_start[i, j])

      elapsed_time_multi_start[i, j] <- system.time({
        fit <- fit_BKP(X, y, m, Xbounds = Xbounds)
      })["elapsed"]
      cat(", BKP_multi_start =", elapsed_time_multi_start[i, j])

      elapsed_time_given_gp[i, j] <- system.time({
        gp <- gp_init(cf = cf_sexp(), lik = lik_binomial())
        gp <- gp_fit(gp, X, y, trials = m)
      })["elapsed"]
      cat(", GP_given =", elapsed_time_given_gp[i, j])

      elapsed_time_optim_gp[i, j] <- system.time({
        gp <- gp_init(cf = cf_sexp(), lik = lik_binomial())
        gp <- gp_optim(gp, X, y, trials = m, verbose = FALSE)
      })["elapsed"]
      cat(", GP_optim =", elapsed_time_optim_gp[i, j], "\n")
    }
  }

  # ---- Compute average time for each n ----
  avg_time <- matrix(NA, nrow = 5, ncol = 5)
  colnames(avg_time) <- paste0("n=", n_vals)
  rownames(avg_time) <- c("given", "single_start", "multi_start", "given_gp", "optim_gp")
  avg_time[1, ] <- colMeans(elapsed_time_given)
  avg_time[2, ] <- colMeans(elapsed_time_single_start)
  avg_time[3, ] <- colMeans(elapsed_time_multi_start)
  avg_time[4, ] <- colMeans(elapsed_time_given_gp)
  avg_time[5, ] <- colMeans(elapsed_time_optim_gp)
}else{
  n_vals <- c(200, 500, 1000, 2000, 5000)
  avg_time <- matrix(c(0.0030, 0.0215, 0.0745, 0.3555, 2.2010, # given
                       0.3095, 1.8790, 6.6080, 24.3170, 144.4830, # single_start
                       5.7440, 31.2160, 132.6125, 547.2550, 3110.1650, # multi_start
                       0.0140, 0.1740, 1.2235, 9.0255, 145.7435, # given_gp
                       0.7500, 8.6575, 73.4760, 471.0000, 6619.0835 # optim_gp
  ), nrow = 5, byrow = TRUE)
  rownames(avg_time) <- c("given", "single_start", "multi_start", "given_gp", "optim_gp")
  colnames(avg_time) <- c("n=200", "n=500", "n=1000", "n=2000", "n=5000")
}

pdf("elapsed_time.pdf", width = 10, height = 5)
par(mfrow = c(1, 2), mar = c(4,4,3,1))  # Adjust margins for better appearance
cols <- c("blue", "skyblue", "darkblue", "red", "orange")
pchs <- c(16, 17, 16, 15, 15)

# Left plot: given and given_gp
plot(n_vals, avg_time["given", ], log = "xy", type = "b", col = cols[1], pch = pchs[1], lwd = 2,
     xlab = "n", ylab = "Time (seconds, log scale)",
     main = "Fixed Hyperparameter",
     ylim = range(avg_time[c("given", "given_gp"), ]))
lines(n_vals, avg_time["given_gp", ], type = "b", col = cols[4], pch = pchs[4], lwd = 2)

# Complexity reference lines
ref_bkp <- avg_time["given", 1] * (n_vals / n_vals[1])^2
ref_gp  <- avg_time["given_gp", 1] * (n_vals / n_vals[1])^3
lines(n_vals, ref_bkp, col = "blue", lty = 2)
lines(n_vals, ref_gp, col = "red", lty = 4)

legend("topleft", bty = "n",
       legend = c("BKP", "LGP", "O(n^2)", "O(n^3)"),
       col = c(cols[c(1,4)], "blue", "red"),
       pch = c(pchs[c(1,4)], NA, NA),
       lty = c(1,1,2,4),
       lwd = c(2,2,1,1),
       bg = "white")

# Right plot: single_start, multi_start, optim_gp
plot(n_vals, avg_time["single_start", ], log = "xy", type = "b", col = cols[2], pch = pchs[2], lwd = 2,
     xlab = "n", ylab = "Time (seconds, log scale)",
     main = "Optimization-based Methods",
     ylim = range(avg_time[c("single_start", "multi_start", "optim_gp"), ]))
lines(n_vals, avg_time["multi_start", ], type = "b", col = cols[3], pch = pchs[3], lwd = 2)
lines(n_vals, avg_time["optim_gp", ], type = "b", col = cols[5], pch = pchs[5], lwd = 2)

# Complexity reference lines
ref_bkp2 <- avg_time["single_start", 1] * (n_vals / n_vals[1])^2
ref_gp2  <- avg_time["optim_gp", 1] * (n_vals / n_vals[1])^3
lines(n_vals, ref_bkp2, col = "skyblue", lty = 2)
lines(n_vals, ref_gp2, col = "orange", lty = 4)

legend("topleft", bty = "n",
       legend = c("BKP (single_start)", "BKP (multi_start)", "LGP", "O(n^2)", "O(n^3)"),
       col = c(cols[c(2,3,5)], "skyblue", "orange"),
       pch = c(pchs[c(2,3,5)], NA, NA),
       lty = c(1,1,1,2,4),
       lwd = c(2,2,2,1,1),
       bg = "white")

# Restore default plotting parameters (optional)
par(mfrow = c(1,1))
dev.off()

#-------------------------- Classification ------------------------
#-------------------------- Example 4 ----------------------------
set.seed(123)
# Data points
n <- 250
n_train <- 200
n_test <- 50
data <- mlbench.spirals(n, cycles = 2, sd = 0.05)

X_train <- data$x[1:n_train, ]
y_train <- as.numeric(data$classes[1:n_train]) - 1 # Convert to 0/1 for BKP
X_test <- data$x[(n_train + 1):n, ]
y_test <- as.numeric(data$classes[(n_train + 1):n]) - 1 # Convert to 0/1 for BKP
m <- rep(1, n_train)
Xbounds <- rbind(c(-1.7, 1.7), c(-1.7, 1.7))

# Fit BKP model
BKP_model_Class <- fit_BKP(
  X_train, y_train, m, Xbounds = Xbounds,
  prior = "fixed", r0 = 0.1, loss = "log_loss")

# Predict via BKP
prediction <- predict(BKP_model_Class, X_test)

# Plot ROC
pdf("ex4roc.pdf", width = 6, height = 6)
roc_curve_BKP <- roc(y_test, prediction$mean)
plot(roc_curve_BKP, main = paste("ROC curve for BKP (AUC = ", round(auc(roc_curve_BKP), 3), ")", sep = ""), col = "#1c61b6")
dev.off()
pdf("ex4.pdf", width = 13, height = 6)
plot(BKP_model_Class)
dev.off()

# Fit LGP model with exponential approx
gp <- gp_init(cf = cf_sexp(), lik = lik_bernoulli())
gp <- gp_optim(gp, X_train, y_train, method = method_full(),
               approx = approx_ep(), verbose = FALSE)

# Predict via LGP
prediction_gp <- gp_pred(gp, as.matrix(X_test), transform = TRUE)

# Plot ROC
pdf("ex4rocLGP.pdf", width = 6, height = 6)
roc_curve_LGP <- roc(y_test, prediction_gp$mean)
plot(roc_curve_LGP, main = paste("ROC curve for LGP (AUC = ", round(auc(roc_curve_LGP), 3), ")", sep = ""), col = "#1c61b6")
dev.off()

# Generate the grid for predictions
x1_seq <- seq(-1.7, 1.7, length.out = 80)
x2_seq <- seq(-1.7, 1.7, length.out = 80)
grid <- expand.grid(x1 = x1_seq, x2 = x2_seq)

# Predict the grid for plot
prediction <- gp_pred(gp, as.matrix(grid), transform = TRUE, var = TRUE)
df <- data.frame(x1 = grid$x1, x2 = grid$x2,
                 Mean = prediction$mean,
                 Variance = prediction$var)
p1 <- BKP:::my_2D_plot_fun("Mean", "Predictive Mean", df, X = X_train, y = y_train)
p2 <- BKP:::my_2D_plot_fun("Variance", "Predictive Variance", df, X = X_train, y = y_train)

# Plot results of LGP
pdf("ex4LGP.pdf", width = 13, height = 6)
gridExtra::grid.arrange(p1, p2, ncol = 2)
dev.off()

#-----------------------------------------------------------------
# ============================================================== #
# ========================= DKP Examples ======================= #
# ============================================================== #
#-------------------------- 1D Example ---------------------------
#-------------------------- Example 5 ----------------------------
set.seed(123)
# Define true class probability function (3-class)
true_pi_fun <- function(X) {
  p1 <- 1/(1+exp(-3*X))
  p2 <- (1 + exp(-X^2) * cos(10 * (1 - exp(-X)) / (1 + exp(-X)))) / 2
  return(matrix(c(p1/2, p2/2, 1 - (p1+p2)/2), nrow = length(p1)))
}
# Data points
n <- 30
Xbounds <- matrix(c(-2, 2), nrow = 1)
X <- lhs(n = n, rect = Xbounds)
true_pi <- true_pi_fun(X)
m <- sample(150, n, replace = TRUE)
Y <- t(sapply(1:n, function(i) rmultinom(1, size = m[i], prob = true_pi[i, ])))

# Fit DKP model
DKP_model_1D <- fit_DKP(X, Y, Xbounds = Xbounds)

# New data points
Xnew <- matrix(seq(-2,2, length = 100), ncol=1)
true_pi <- true_pi_fun(Xnew)

# Plot results
pdf("ex5.pdf", width = 8,height = 8)
plot(DKP_model_1D)
plot(Xnew, true_pi[, 1], type = "l", col = "black",
     xlab = "x", ylab = "Probability", ylim = c(0, 1),
     main = "True Probability", lwd = 2)
lines(Xnew, true_pi[, 2], col = "red", lwd = 2)
lines(Xnew, true_pi[, 3], col = "blue", lwd = 2)
legend("topright",
       bty = "n",
       legend = c("Class 1", "Class 2", "Class 3"),
       col = c("black", "red", "blue"), lty = 1, lwd = 2)
dev.off()

#-------------------------- 2D Example ---------------------------
#-------------------------- Example 6 ----------------------------
set.seed(123)
# Define 2D latent function and probability transformation (3-class)
true_pi_fun <- function(X){
  if(is.null(nrow(X))) X <- matrix(X, nrow=1)
  m <- 8.6928
  s <- 2.4269
  x1 <- 4*X[,1]- 2
  x2 <- 4*X[,2]- 2
  a <- 1 + (x1 + x2 + 1)^2 *
    (19- 14*x1 + 3*x1^2- 14*x2 + 6*x1*x2 + 3*x2^2)
  b <- 30 + (2*x1- 3*x2)^2 *
    (18- 32*x1 + 12*x1^2 + 48*x2- 36*x1*x2 + 27*x2^2)
  f <- (log(a*b)- m)/s
  p1 <- pnorm(f)  # Transform to probability
  p2 <- sin(pi * X[,1]) * sin(pi * X[,2])
  return(matrix(c(p1/2, p2/2, 1 - (p1+p2)/2), nrow = length(p1)))
}
# Data points
n <- 100
Xbounds <- matrix(c(0, 0, 1, 1), nrow = 2)
X <- lhs(n = n, rect = Xbounds)
true_pi <- true_pi_fun(X)
m <- sample(150, n, replace = TRUE)
Y <- t(sapply(1:n, function(i) rmultinom(1, size = m[i], prob = true_pi[i, ])))

# Fit DKP model
DKP_model_2D <- fit_DKP(X, Y, Xbounds=Xbounds)

# Plot results
pdf(file = "ex6_class%d.pdf", width = 9, height = 8, onefile = FALSE)
plot(DKP_model_2D)
dev.off()

# New data points
Xnew1 <- seq(Xbounds[1,1], Xbounds[1,2], length.out = 200)
Xnew2 <- seq(Xbounds[2,1], Xbounds[2,2], length.out = 200)
Xnew <- expand.grid(Xnew1 = Xnew1, Xnew2 = Xnew2)
true_pi <- true_pi_fun(Xnew)
df <- data.frame(x1 = Xnew$Xnew1, x2 = Xnew$Xnew2,
                 True1 = true_pi[,1],
                 True2 = true_pi[,2],
                 True3 = true_pi[,3])
# Plot True distribution
pdf("ex6_class1_true.pdf", width = 4.5, height = 4)
print(BKP:::my_2D_plot_fun("True1", title = "True Probability", data = df))
dev.off()

pdf("ex6_class2_true.pdf", width = 4.5, height = 4)
print(BKP:::my_2D_plot_fun("True2", title = "True Probability", data = df))
dev.off()

pdf("ex6_class3_true.pdf", width = 4.5, height = 4)
print(BKP:::my_2D_plot_fun("True3", title = "True Probability", data = df))
dev.off()


#-------------------------- Classification ------------------------
#-------------------------- Example 7 ----------------------------
set.seed(123)
data(iris)
X <- as.matrix(iris[, 1:2])
Xbounds <- rbind(c(4.2, 8), c(1.9, 4.5))
labels <- iris$Species
Y <- model.matrix(~ labels - 1) # expand factors to a set of dummy variables

train_indices <- sample(1:nrow(iris), 0.7 * nrow(iris))
X_train <- X[train_indices, ]
Y_train <- Y[train_indices, ]
labels_train <- labels[train_indices]

X_test <- X[-train_indices, ]
Y_test <- Y[-train_indices, ]
labels_test <- labels[-train_indices]

# Fit DKP model
DKP_model_Class <- fit_DKP(
  X_train, Y_train, Xbounds = Xbounds, loss = "log_loss",
  prior = "fixed", r0 = 0.1, p0 = rep(1/3, 3))
pdf("ex7.pdf", width = 13, height = 6)
plot(DKP_model_Class)
dev.off()

# Test
dkp_pred_probs <- predict(DKP_model_Class, X_test)$mean
class_levels <- levels(labels_test)
colnames(dkp_pred_probs) <- class_levels
multiclass_roc_dkp <- multiclass.roc(labels_test, dkp_pred_probs)

all_rocs <- list()
for (class_name in class_levels) {
  # convert to binary label
  labels_binary <- ifelse(labels_test == class_name, 1, 0)
  probabilities <- dkp_pred_probs[, class_name]
  roc_curve <- roc(labels_binary, probabilities)
  all_rocs[[class_name]] <- roc_curve
}

# Plot ROC curve
pdf("ex7roc.pdf", width = 6, height = 6)
plot(all_rocs[[1]], col = "blue", lwd = 2, lty = 1,
     main = paste("One-vs-Rest ROC curve for BKP (AUC =", round(auc(multiclass_roc_dkp), 3), ")", sep = ""))
lines(all_rocs[[2]], col = "red", lwd = 2, lty = 2)
lines(all_rocs[[3]], col = "black", lwd = 2, lty = 4)
legend("bottomright",
       legend = class_levels,
       col = c("blue", "red", "black"),
       lwd = 2,
       lty = c(1,2,4),
       cex = 1.2)
dev.off()

### LGP model
iris_data <- data.frame(
  Sepal.Length = iris$Sepal.Length,
  Sepal.Width = iris$Sepal.Width,
  Species = iris$Species
)

iris_train <- iris_data[train_indices, ]
iris_test <- iris_data[-train_indices, ]

# Fit LGP model
gausspr_model <- gausspr(Species ~ ., data = iris_train,
                         kernel = "rbfdot", kpar = "automatic")

# Test LGP
lgp_pred_probs <- predict(gausspr_model, newdata = iris_test,
                          type = "probabilities")
multiclass_roc_lgp <- multiclass.roc(iris_test$Species, lgp_pred_probs)

# Plot ROC
all_rocs <- list()
for (class_name in class_levels) {
  # binary label
  labels_binary <- ifelse(labels_test == class_name, 1, 0)
  probabilities <- lgp_pred_probs[, class_name]
  roc_curve <- roc(labels_binary, probabilities)
  all_rocs[[class_name]] <- roc_curve
}

pdf("ex7rocLGP.pdf", width = 6, height = 6)
plot(all_rocs[[1]], col = "blue", lwd = 2, lty = 1,
     main = paste("One-vs-Rest ROC curve for LGP (AUC =", round(auc(multiclass_roc_lgp), 3), ")", sep = ""))
lines(all_rocs[[2]], col = "red", lwd = 2, lty = 2)
lines(all_rocs[[3]], col = "black", lwd = 2, lty = 4)
legend("bottomright",
       legend = class_levels,
       col = c("blue", "red", "black"),
       lty = c(1,2,4),
       lwd = 2,
       cex = 1.2)
dev.off()

# New data points for plot
grid <- expand.grid(
  Sepal.Length = seq(Xbounds[1, 1], Xbounds[1, 2], length.out = 80),
  Sepal.Width = seq(Xbounds[2, 1], Xbounds[2, 2], length.out = 80)
)
grid_predictions <- predict(gausspr_model, newdata = grid, type = "prob")
class <- max.col(grid_predictions)
df<- data.frame(x1 = grid$Sepal.Length,x2=grid$Sepal.Width,
                class = class,
                max_prob = apply(grid_predictions, 1, max))
p1 <- BKP:::my_2D_plot_fun_class("class", "Predicted Classes", df, X_train, Y_train)
p2 <- BKP:::my_2D_plot_fun_class("max_prob", "Maximum Predicted Probability", df, X_train, Y_train, classification = FALSE)
pdf("ex7LGP.pdf", width = 13, height = 6)
grid.arrange(p1, p2, ncol = 2)
dev.off()
```

## BKP Model  {#sec:example_bkp}

This subsection presents detailed illustrative examples based on the BKP model. Examples for the DKP model are provided in Section \ref{sec:example_dkp}.

\paragraph{Example 1} Let $x\in[-2,2]$, and suppose the true Bernoulli probability function is given by
\begin{equation}\label{eq:bkp-1d-generate-function-1}
	\pi_{1}(x) = \frac{1}{1+e^{-3x}},
\end{equation}
which is referred as the function \code{true_pi_fun} in the code below.
We aim to fit the BKP model based on seven input locations that are uniformly distributed over $[-2,2]$, with each location associated with a binomial observation having a maximum trial count of 100. The input locations are generated using the \fct{lhs} function from the \proglang{R} package \pkg{tgp} \citep{Gramacy2010tgp2}. The following \proglang{R} code illustrates how to simulate the data and fit the BKP model using the \fct{fit\_BKP} function.

\begin{CodeChunk}
\begin{CodeInput}    
R> n <- 7
R> Xbounds <- matrix(c(-2,2), nrow = 1) 
R> X <- lhs(n = n, rect = Xbounds) 
R> true_pi <- true_pi_fun(X) 
R> m <- sample(100, n, replace = TRUE)
R> y <- rbinom(n, size = m, prob = true_pi) 
R> BKP_model_1D_1 <- fit_BKP(X, y, m, Xbounds = Xbounds)
\end{CodeInput}
\end{CodeChunk}

The estimates of the parameters of the fitted BKP model can be displayed using the \fct{print} function:

\begin{CodeChunk}
\begin{CodeInput} 
R> print(BKP_model_1D_1)
\end{CodeInput}
\begin{CodeOutput}  

       Beta Kernel Process (BKP) Model    

Number of observations (n):  7
Input dimensionality (d):    1
Kernel type:                 gaussian
Optimized kernel parameters: 0.1748
Minimum achieved loss:       0.01165
Loss function:               brier
Prior type:                  noninformative
\end{CodeOutput}
\end{CodeChunk}

The \code{BKP_model_1D_1} object can be used for visualization and prediction over a grid of input values via the \fct{plot} and \fct{simulate} methods:

\begin{CodeChunk}
\begin{CodeInput}  
R> plot(BKP_model_1D_1)
R> Xnew = matrix(seq(-2, 2, length = 100), ncol = 1)
R> sim <- simulate(BKP_model_1D_1, nsim = 3, Xnew = Xnew)
\end{CodeInput} 
\end{CodeChunk}

\begin{figure}[!t]
	\centering
	\begin{subfigure}{0.49\textwidth}
		\includegraphics[width=\linewidth]{ex1.pdf}
		\caption{Posterior mean and 95\% CI of Example 1}
		\label{fig:ex1}
	\end{subfigure}
	\hfill 
	\begin{subfigure}{0.49\textwidth}
		\includegraphics[width=\linewidth]{ex1_sim.pdf}
		\caption{Posterior samples of Example 1}
		\label{fig:ex1_sim}
	\end{subfigure}
	\medskip
	\begin{subfigure}{0.49\textwidth}
		\includegraphics[width=\linewidth]{ex2.pdf}
		\caption{Posterior mean and 95\% CI of Example 2}
		\label{fig:ex2}
	\end{subfigure}
	\hfill 
	\begin{subfigure}{0.49\textwidth}
		\includegraphics[width=\linewidth]{ex2_sim.pdf}
		\caption{Posterior samples of Example 2}
		\label{fig:ex2_sim}
	\end{subfigure} 
	
	\caption{Posterior inference and simulation results from the fitted BKP models}
	\label{fig:BKP_1D}
\end{figure}

Figure \ref{fig:BKP_1D} presents two views of the model output. Panel~(a) shows the posterior mean estimate of the probability function $\pi(x)$ (blue), along with a 95\% credible interval (gray band), observed proportions (red dots), and the true underlying probability function (black). Panel~(b) presents three posterior sample curves generated using the \fct{simulate} method, illustrating the variability of the estimated probability surface.

We continue with the same probability function of this example to show the impact of the global precision parameters, $r_0$, in classification tasks. However, the sample size is changed from 7 to 20 for classification. The following \proglang{R} code generates the label of response and fit the BKP model for classification with $r_0$ being 0.01 and 2. 

\begin{CodeChunk}
\begin{CodeInput}  
R> # Fit BKP model with r0 = 0.01
R> BKP_model_1D_1_class_1 <- fit_BKP(
+    X, y, m, Xbounds = Xbounds,
+    prior = "fixed", r0 = 0.01, loss = "log_loss")
R> # Fit BKP model with r0 = 2
R> BKP_model_1D_1_class_2 <- fit_BKP(
+    X, y, m, Xbounds = Xbounds,
+    prior = "fixed", r0 = 2, loss = "log_loss")
\end{CodeInput}
\end{CodeChunk} 
 
Figure \ref{fig:ex1_class_001} shows a sigmoidal curve with steep slope around zero (in solid line) and a narrow 95\% credible interval (in grey band) under $r_0 = 0.01$, indicating a decisive and confident classification boundary. In contrast, Figure \ref{fig:ex1_class_2} displays a sine-shape curve with a much wider credible interval, reflecting greater uncertainty and less effective separation. This demonstrates the preference of a small value of $r_0$ value for classification task. 
 
\begin{figure}[!t]
	\centering
	\begin{subfigure}{0.49\textwidth}
		\includegraphics[width=\linewidth]{ex1class001.pdf}
		\caption{$r_0=0.01$}
		\label{fig:ex1_class_001}
	\end{subfigure}
	\hfill 
	\begin{subfigure}{0.49\textwidth}
		\includegraphics[width=\linewidth]{ex1class2.pdf}
		\caption{$r_0=2$}
		\label{fig:ex1_class_2}
	\end{subfigure}
	\caption{Posterior mean and 95\% CI of Example 1 (classification task) under $r_0=0.01$ and 2}
	\label{fig:BKP_1D_class}
\end{figure}

\paragraph{Example 2} The first example is essentially a generalized linear model with a smooth logit link, and thus poses limited modeling complexity. To demonstrate the capability of the BKP model in handling more challenging classification structures, we consider a second example with a highly nonlinear underlying probability surface. Define the true Bernoulli probability as
\begin{equation}\label{eq:bkp-1d-generate-function-2}
	\pi_{2}(x) = \frac{1}{2}\left[ 1+e^{-x^{2}}\cos \left( 10\frac{1-e^{-x}}{1+e^{-x}} \right) \right],
\end{equation}
where $x \in [-2,2]$  \citep{Goetschalckx:2011}. This example involves rapid local oscillations and strong nonlinearity, making it substantially more difficult to fit than Example~1. Here, we increase the number of locations to 30.

\begin{CodeChunk}
\begin{CodeInput} 
R> n <- 30
R> Xbounds <- matrix(c(-2,2), nrow = 1)
R> X <- lhs(n = n, rect = Xbounds)
R> true_pi <- true_pi_fun(X)
R> m <- sample(100, n, replace = TRUE)
R> y <- rbinom(n, size = m, prob = true_pi) 
R> BKP_model_1D_2 <- fit_BKP(X, y, m, Xbounds = Xbounds)
\end{CodeInput}
\end{CodeChunk}

The results are shown in panels~(c) and (d) of Figure \ref{fig:BKP_1D}. Panel~(c) demonstrates that the BKP model accurately recovers the highly nonlinear pattern, while panel~(d) illustrates posterior variability through three representative realizations of $\pi(x)$.

\paragraph{Example 3} We now consider a two-dimensional test function to further illustrate the modeling capabilities of the \pkg{BKP} package. Let $\vec{x} \in [0,1]^2$, and define the latent surface using a re-scaled version of the Goldstein–Price function \citep{Picheny2013Benchmark}:
\begin{align*}
	f(\vec{x})=& \frac{\log[\{1+a(\vec{x})\}\{30+b(\vec{x})\}-8.6928]}{2.4269}, \quad \text{with}\\
	a(\vec{x}) =& \left(4 x_1 + 4 x_2 - 3 \right)^2  \times  \notag  \\
	& \{75 - 56 \left(x_1 + x_2 \right) + 3\left(4 x_1 - 2 \right)^2 + 6\left(4 x_1 - 2 \right)\left(4 x_2 - 2 \right) + 3\left(4 x_2 - 2 \right)^2\}, \notag \\
	b(\vec{x}) =& \left(8 x_1 - 12 x_2 +2 \right)^2 \times \notag \\
	& \{-14 - 128 x_1 + 12\left(4 x_1 - 2 \right)^2 + 192 x_2 - 36\left(4 x_1 - 2 \right)\left(4 x_2 - 2 \right) + 27\left(4 x_2 - 2 \right)^2 \}.\notag 
\end{align*} 

The true Bernoulli probability surface is then defined by 
\begin{equation}\label{eq:Goldstein-Price}
	\pi_3(\vec{x}) = \Phi\{f(\vec{x})\}, 
\end{equation}
where $\Phi(\cdot)$ is the cumulative distribution function of the standard normal distribution. This formulation produces a smooth yet highly non-linear response surface, providing a challenging test scenario for probabilistic modeling.

To construct the training data, we generate a LHD of size $100$ over $[0,1]^2$. Each location is associated with a binomial observation whose number of trials is randomly drawn from $\{1, \dots, 100\}$. The following \proglang{R} code demonstrates the data simulation and model fitting using the \fct{fit\_BKP} function:

\begin{CodeChunk}
\begin{CodeInput}
R> n <- 100
R> Xbounds <- matrix(c(0, 0, 1, 1), nrow = 2)
R> X <- lhs(n = n, rect = Xbounds)
R> true_pi <- true_pi_fun(X)
R> m <- sample(100, n, replace = TRUE)
R> y <- rbinom(n, size = m, prob = true_pi) 
R> BKP_model_2D <- fit_BKP(X, y, m, Xbounds=Xbounds) 
R> print(BKP_model_2D)
\end{CodeInput}
\begin{CodeOutput} 

       Beta Kernel Process (BKP) Model    

Number of observations (n):  100
Input dimensionality (d):    2
Kernel type:                 gaussian
Optimized kernel parameters: 0.1112, 0.0680
Minimum achieved loss:       0.01041
Loss function:               brier
Prior type:                  noninformative
\end{CodeOutput}
\end{CodeChunk} 

\begin{figure}[!t]
	\centering 
	\begin{subfigure}{0.45\textwidth}
		\includegraphics[width=\linewidth]{ex3_true.pdf} 
	\end{subfigure}   
	\medskip  
	\begin{subfigure}{\textwidth}
		\includegraphics[width=\linewidth]{ex3.pdf}  
	\end{subfigure} 
	\caption{Comparison of the true probability surface and BKP-based posterior summaries of Example~3} 
	\label{fig:BKP_2D}
\end{figure}

Figure \ref{fig:BKP_2D} is obtained by the following code, which is the heatmaps with contour lines for the true probability and the estimated. 

\begin{CodeChunk}
\begin{CodeInput}  
R> plot(BKP_model_2D) 
\end{CodeInput} 
\end{CodeChunk}

The results of model prediction are as follows:

\begin{CodeChunk}
\begin{CodeInput}
R> Xnew <- lhs(n = 10, rect = Xbounds)
R> predict(BKP_model_2D, Xnew)
\end{CodeInput}
\begin{CodeOutput}
Prediction results on new data (Xnew).
Total number of prediction points: 10 

Preview of predictions for new data (first 6 of 10 points):
     x1     x2   mean variance 2.5% quantile 97.5% quantile
 0.8249 0.0605 0.7106   0.0012        0.6394         0.7770
 0.3784 0.5314 0.3324   0.0009        0.2754         0.3919
 0.0566 0.1437 0.5519   0.0032        0.4409         0.6603
 0.5402 0.8851 0.7994   0.0033        0.6756         0.8993
 0.7031 0.7113 0.3282   0.0009        0.2707         0.3883
 0.2709 0.6661 0.6178   0.0013        0.5460         0.6871
 ...
\end{CodeOutput}
\end{CodeChunk} 

We continue with Example 3 to demonstrate the scalability of the BKP model in comparison to the logistic Gaussian process (LGP) model \citep{Rasmussen2006GPML}. Let the sample size range from 200 to 5000. Two scenarios for hyperparameter determination are considered, namely, i) one scenario with fixed value of \code{theta = 1}, and ii) the other scenario with hyperparameters optimized (through \code{n_multi_start = 1} or \code{n_multi_start = NULL}). 

The LGP model was implemented using the \fct{gp\_fit} and \fct{gp\_optim} functions from the \pkg{gplite} package \citep{Piironen2022gplite}. Notably, the \pkg{gplite} package employs a multi-start optimization strategy that is only triggered upon optimization failure; therefore, the additional computational cost of multiple restarts was not included in our timing measurements.

All experiments were conducted on a workstation equipped with an Intel(R) Xeon(R) W-2235 CPU @ 3.80GHz (12 cores) and 16 GB RAM. The computation times reported here are averages over 20 independent repetitions to mitigate the effects of variability due to random initialization and computational fluctuations.

The average computation times for both models under these scenarios are presented in Figure \ref{fig:elapsed_time}. The observed computational costs align well with the theoretical complexity predictions, confirming the expected $\mathcal{O}(n^2)$ scalability for BKP and $\mathcal{O}(n^3)$ for the LGP model. Additionally, the computational cost of the multi-start optimization approach for BKP is roughly proportional to the number of restarts, which was set to $10d$. For the current setting, this corresponds to approximately 20 times the cost of the single-start method, reflecting the increased complexity incurred by multiple initializations. This trade-off between computational expense and potentially improved optimization robustness should be carefully considered when selecting an appropriate strategy.
 


\begin{figure}[!t]
	\centering 
	\includegraphics[width=\linewidth]{elapsed_time.pdf}  
	\caption{Comparison of computation times (in log scale) between BKP and LGP methods: (a) fixed hyperparameter; (b) optimization-based methods}
	\label{fig:elapsed_time}
\end{figure}

\paragraph{Example 4}  
We next consider a binary classification task using the \emph{Two Spirals} dataset \citep{Chalup2007spirals}, a well-known benchmark consisting of two intertwined spirals in a bounded two-dimensional input space. This dataset is particularly challenging due to the complex, non-linearly separable class structure.

We generate $n = 250$ observations using the \fct{mlbench.spirals} function from the \proglang{R} package \pkg{mlbench} \citep{Leisch2024mlbench}, with two complete rotations and additive Gaussian noise of standard deviation \code{sd = 0.05}. The inputs $\vec{x}$ are constrained to the domain $[-1.7, 1.7]^2$, and the binary class labels are encoded as 0 and 1. We fit the BKP model using a fixed prior specification with \code{r_0 = 0.1} and \code{p_0 = 0.5}.

\begin{CodeChunk}
\begin{CodeInput}
R> n <- 250
R> n_train <- 200
R> n_test <- 50
R> data <- mlbench.spirals(n, cycles = 2, sd = 0.05)
R> X_train <- data$x[1:n_train, ]
R> y_train <- as.numeric(data$classes[1:n_train]) - 1 # Convert to 0/1 for BKP
R> X_test <- data$x[(n_train + 1):n, ]
R> y_test <- as.numeric(data$classes[(n_train + 1):n]) - 1  
R> m <- rep(1, n_train)
R> Xbounds <- rbind(c(-1.7, 1.7), c(-1.7, 1.7)) 
R> BKP_model_Class <- fit_BKP(
+    X_train, y_train, m, Xbounds = Xbounds,
+    prior = "fixed", r0 = 0.1, loss = "log_loss") 
R> prediction <- predict(BKP_model_Class, X_test)
\end{CodeInput}
\end{CodeChunk} 

\begin{figure}[!t]
	\centering 
	\begin{subfigure}{\textwidth}
		\includegraphics[width=\linewidth]{ex4.pdf}  
		\caption{BKP}
	\end{subfigure} 
	\medspace
	\begin{subfigure}{\textwidth}
		\includegraphics[width=\linewidth]{ex4LGP.pdf}  
		\caption{LGP}
	\end{subfigure}    
  	\caption{Predictions of the two spirals by BKP and LGP models in terms of posterior predictive mean (left panel) and predictive variance (right panel) of class probabilities, where color shading represents the estimated success probability, circles and crosses indicate training observations from the two classes}
	\label{fig:BKP_spiral}
\end{figure}

Figure \ref{fig:BKP_spiral} presents the predictive mean and variance surfaces of the class probabilities. The upper left panel shows that the BKP model accurately captures the intricate spiral structure, with smoothly varying predicted probabilities that delineate the nonlinear decision boundary. The upper right panel displays the associated predictive uncertainty, which is highest near the boundary regions between the spirals. Circles and crosses indicate training observations from the two respective classes. These results illustrate the capacity of BKP to flexibly model complex classification boundaries while providing coherent uncertainty quantification. 

We compare the prediction performance with the LGP model implemented by R package \pkg{gplite}. 
\begin{CodeChunk}
\begin{CodeInput}
R> gp <- gp_init(cf = cf_sexp(), lik = lik_bernoulli())
R> gp <- gp_optim(gp, X_train, y_train, method = method_full(),
+                 approx = approx_ep(), verbose = FALSE) 
R> prediction_gp <- gp_pred(gp, as.matrix(X_test), transform = TRUE)
\end{CodeInput}
\end{CodeChunk} 

The lower panel of Figure \ref{fig:BKP_spiral} indicates that the LGP model encounters difficulty in capturing the intricate geometry of the Two Spirals dataset. Although it roughly follows the spiral structure, its predictive mean exhibits a more irregular and less sharply defined decision boundary compared to the BKP model. This pattern suggests that the LGP model may be more prone to overfitting or less capable of generalizing effectively. The corresponding predictive variance plots reinforce this interpretation when viewed separately in regions with and without training data. In regions where data are available, the LGP model exhibits elevated variance (brighter regions) across a broader portion of the input space indicating greater uncertainty in areas critical for classification. In contrast, the BKP model maintains lower variance concentrated more closely around the decision boundary, reflecting higher confident predictions. Interestingly, in regions without training data, the LGP model’s uncertainty is not necessarily the highest; instead, it sometimes shows comparatively lower variance than the BKP model. This suggests that the LGP model may over-smooth or fail to express appropriate uncertainty far from observed data, whereas BKP better captures uncertainty behavior in unobserved regions. The receiver operating characteristic (ROC) curves and the corresponding area under the curve (AUC) values in Figure \ref{fig:ex4-roc} further corroborate this conclusion: BKP achieves an AUC of 0.926, compared to 0.889 for LGP, reflecting superior discriminatory performance and a better balance between sensitivity and specificity. 

\begin{figure}[!t]
	\centering 
	\begin{subfigure}{0.48\textwidth}
		\includegraphics[width=\linewidth]{ex4roc.pdf}  
		\caption{BKP}
	\end{subfigure} 
	\begin{subfigure}{0.48\textwidth}
		\includegraphics[width=\linewidth]{ex4rocLGP.pdf}  
		\caption{LGP}
	\end{subfigure}    
  	\caption{ROC curves and their respective AUCs for BKP and LGP models}
	\label{fig:ex4-roc}
\end{figure}

## DKP Model {#sec:example_dkp}
 
\paragraph{Example 5}
Consider a one-dimensional three-class classification problem. The input is defined on the interval $x \in [-2, 2]$, and the true class probability vector is given by
\begin{equation*} 
	\vec{\pi}(x) = \left[\frac{\pi_{1}(x)}{2}, \frac{\pi_{2}(x)}{2},1-\frac{\pi_{1}(x)}{2}-\frac{\pi_{2}(x)}{2}\right]^\top,
\end{equation*}
where $\pi_{1}(x)$ and $\pi_{2}(x)$ are smooth functions defined in \eqref{eq:bkp-1d-generate-function-1} and \eqref{eq:bkp-1d-generate-function-2}, respectively.

We generate $n = 30$ input locations using Latin hypercube sampling over the interval $[-2, 2]$. At each location, the response is a multinomial vector with probability $\vec{\pi}(x)$ and a random total count sampled from $\{1, \dots, 150\}$. The DKP model is then fitted using the \fct{fit\_DKP} function and visualized with the \fct{plot} method:

\begin{CodeChunk}
\begin{CodeInput}  
R> n <- 30
R> Xbounds <- matrix(c(-2, 2), nrow = 1)
R> X <- lhs(n = n, rect = Xbounds)
R> true_pi <- true_pi_fun(X)
R> m <- sample(150, n, replace = TRUE)
R> Y <- t(sapply(1:n, function(i) rmultinom(1, size = m[i], prob = true_pi[i, ]))) 
R> DKP_model_1D <- fit_DKP(X, Y, Xbounds = Xbounds)
R> plot(DKP_model_1D) 
\end{CodeInput}
\end{CodeChunk} 

\begin{figure}[!t]
	\centering 
	\includegraphics[width=\linewidth]{ex5.pdf}  
	\caption{Posterior inference from the fitted DKP model for a one-dimensional three-class problem} 
	\label{fig:ex5}
\end{figure} 

Figure \ref{fig:ex5} shows that the DKP model accurately recovers the true class probability functions and provides meaningful uncertainty quantification. The predictive mean curves align with the true underlying structure, and the shaded bands reflect posterior uncertainty.

\paragraph{Example 6} Consider a two-dimensional three-class classification problem. Let $\vec{x} = [x_{1}, x_{2}]^\top \in [0,1]^2$, and define the true class probability function as
\begin{align*} 
	\vec{\pi}(\vec{x})= \left[\frac{\pi_{3}(\vec{x})}{2}, \; \frac{\pi_{4}(\vec{x})}{2}, \; 1-\frac{\pi_{3}(\vec{x})}{2}- \frac{\pi_{4}(\vec{x})}{2}\right]^\top,
\end{align*} 
where  $\pi_{3}(\vec{x})$ is defined in \eqref{eq:Goldstein-Price}, and
\[
	\pi_{4}(\vec{x}) = \sin(\pi x_1) \sin(\pi x_2).
\]

We generate $n = 100$ input points using Latin hypercube sampling over $[0,1]^2$. 
At each location, the response is a multinomial vector with probability $\vec{\pi}(\vec{x})$ and a total count randomly sampled from $\{1, \dots, 150\}$.
We fit the DKP model using \fct{fit\_DKP} and visualize the results with the \fct{plot} method:

\begin{CodeChunk}
\begin{CodeInput} 
R> n <- 100
R> Xbounds <- matrix(c(0, 0, 1, 1), nrow = 2)
R> X <- lhs(n = n, rect = Xbounds)
R> true_pi <- true_pi_fun(X)
R> m <- sample(150, n, replace = TRUE)
R> Y <- t(sapply(1:n, function(i) rmultinom(1, size = m[i], prob = true_pi[i, ]))) 
R> DKP_model_2D <- fit_DKP(X, Y, Xbounds=Xbounds) 
\end{CodeInput}
\end{CodeChunk} 

\begin{figure}[!t]
	\centering 
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\linewidth]{ex6_class1_true.pdf} 
	\end{subfigure}   
	\medspace 
	\begin{subfigure}{\textwidth}
		\includegraphics[width=\linewidth]{ex6_class1.pdf}  
	\end{subfigure} 
	\caption{Class 1: true distribution (top) and DKP model posterior summaries (bottom)} 
	\label{fig:ex6-class1}
\end{figure} 

\begin{figure}[!t]
	\centering 
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\linewidth]{ex6_class2_true.pdf} 
	\end{subfigure}   
	\medspace   
	\begin{subfigure}{\textwidth}
		\includegraphics[width=\linewidth]{ex6_class2.pdf}  
	\end{subfigure} 
	\caption{Class 2: true distribution (top) and DKP model posterior summaries (bottom)} 
	\label{fig:ex6-class2}
\end{figure} 

\begin{figure}[!t]
	\centering 
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\linewidth]{ex6_class3_true.pdf} 
	\end{subfigure}   
	\medspace 
	\begin{subfigure}{\textwidth}
		\includegraphics[width=\linewidth]{ex6_class3.pdf}  
	\end{subfigure} 
	\caption{Class 3: true distribution (top) and DKP model posterior summaries (bottom)} 
	\label{fig:ex6-class3}
\end{figure}

Figures~\ref{fig:ex6-class1}--\ref{fig:ex6-class3} display the ground truth surfaces and corresponding posterior inference for each of the three classes. In each figure, the top panel shows the true class probability surface, while the bottom row presents the posterior mean, variance, and uncertainty bounds produced by the DKP model. 

\paragraph{Example 7} 
Consider another three-class classification task based on the well-known \emph{Iris} dataset available in \proglang{R}. This classic benchmark dataset contains measurements of three iris species---\emph{setosa}, \emph{versicolor}, and \emph{virginica}---each described by four features: sepal length, sepal width, petal length, and petal width. Due to the overlap in feature space, particularly between \emph{versicolor} and \emph{virginica}, the class boundaries are not linearly separable, making this dataset a standard testbed for evaluating multi-class classification algorithms. 
For visualization purposes, we restrict our analysis to the first two features: sepal length and sepal width. We fit the DKP model using a fixed prior specification with \code{r0 = 0.1} and \code{p0 = rep(1/3, 3)}, assuming equal prior probability for each class.

\begin{CodeChunk}
\begin{CodeInput}
R> data(iris)
R> X <- as.matrix(iris[, 1:2])
R> Xbounds <- rbind(c(4.2, 8), c(1.9, 4.5))
R> labels <- iris$Species
R> Y <- model.matrix(~ labels - 1)  # expand factors to a set of dummy variables
R> train_indices <- sample(1:nrow(iris), 0.7 * nrow(iris))
R> X_train <- X[train_indices, ]
R> Y_train <- Y[train_indices, ]
R> DKP_model_Class <- fit_DKP(
+    X_train, Y_train, Xbounds = Xbounds, loss = "log_loss",
+    prior = "fixed", r0 = 0.1, p0 = rep(1/3, 3))
R> X_test <- X[-train_indices, ]
R> Y_test <- Y[-train_indices, ]
R> dkp_pred_probs <- predict(DKP_model_Class, X_test)$mean
\end{CodeInput}
\end{CodeChunk} 

\begin{figure}[!t]
	\centering 
	\begin{subfigure}{\textwidth}
		\includegraphics[width=\linewidth]{ex7.pdf} 
		\caption{DKP }
	\end{subfigure} 
	\begin{subfigure}{\textwidth}
		\includegraphics[width=\linewidth]{ex7LGP.pdf} 
		\caption{LGP }
	\end{subfigure}  
  	\caption{Classification and uncertainty visualization after applying the DKP model on the \emph{Iris} dataset. The left panel shows the predicted classification regions based on the MAP decision rule, where each color corresponds to one of the three iris species. The training data points are overlaid using shape-coded markers. The right panel visualizes the model's predictive uncertainty using the maximum predicted class probability. Lower values (light-colored regions) indicate higher classification uncertainty, which usually occurs near the class boundaries.}
	\label{fig:ex7}
\end{figure}

Figure \ref{fig:ex7} provides a visual summary of the fitted DKP model. 
The upper left panel displays the MAP classification boundaries, which clearly separate \emph{setosa} from the other two species. In contrast, the boundary between \emph{versicolor} and \emph{virginica} is more intricate, reflecting the known overlap between these two species in sepal measurements.

To further illustrate the model's prediction uncertainty, the upper right panel shows a contour map of the maximum predicted probability $\max_j \pi_j(\vec{x})$ across the input space. 
This diagnostic highlights regions where the classifier is most confident (values near 1) versus uncertain (values close to 1/3), which effectively identifies decision boundaries and ambiguous areas. 
Such visualization aids in understanding the reliability of classification decisions and the structure of the learned decision surfaces.

Again, we compare the performance with LGP model, which is implemented by R package \pkg{kernlab} \citep{Karatzoglou2004kernlab}. We adopt the one-vs-rest (OvR) approach for comparison, which works by training a separate binary classifier for each class. Each classifier is tasked with distinguishing its assigned class from all other classes combined.

\begin{CodeChunk}
\begin{CodeInput}
R> iris_data <- data.frame(
+    Sepal.Length = iris$Sepal.Length,
+    Sepal.Width = iris$Sepal.Width,
+    Species = iris$Species
+  )
R> iris_train <- iris_data[train_indices, ]
R> iris_test <- iris_data[-train_indices, ] 
R> gausspr_model <- gausspr(Species ~ ., data = iris_train, 
+                           kernel = "rbfdot", kpar = "automatic") 
R> lgp_pred_probs <- predict(gausspr_model, newdata = iris_test, 
+                            type = "probabilities")
\end{CodeInput}
\end{CodeChunk} 

It is seen that the decision boundaries in the ``Predicted Classes'' plot (lower left panel) appear less smooth, particularly in the regions between the two clusters of `triangle' and `plus' data points. Additionally, the ``Maximum Predicted Probability'' plot (lower right panel) shows a more fragmented and less uniform confidence landscape. There are pockets of high confidence, but also areas of lower confidence (yellows and greens) scattered over the region, even within what appear to be the core regions of the classes. This indicates that DKP model is more accurate in prediction and more consistent in uncertainty quantification than the LGP model. The corresponding multiclass ROC curves with macro-average AUC (0.936 vs 0.927) are presented in Figure \ref{fig:ex7-roc}.  

\begin{figure}[!t]
	\centering 
	\begin{subfigure}{0.48\textwidth}
		\includegraphics[width=\linewidth]{ex7roc.pdf}  
		\caption{BKP}
	\end{subfigure} 
	\begin{subfigure}{0.48\textwidth}
		\includegraphics[width=\linewidth]{ex7rocLGP.pdf}  
		\caption{LGP}
	\end{subfigure}    
	\caption{One-vs-Rest ROC curves and their respective macro-average AUCs for BKP and LGP models}
	\label{fig:ex7-roc}
\end{figure}


# Working Example: Modeling Loa loa Infection {#sec:Loaloa}

```{r, Working Example, include=FALSE}
# Load the data
data("Loaloa")
# Extract input variables (X), response variable (y), and trial counts (m)
X <- as.matrix(Loaloa[, 1:2])
rownames(X) <- NULL
y <- Loaloa$npos
m <- Loaloa$ntot

# Randomly split into training (70%) and testing (30%) sets
set.seed(123)
train_idx <- sample(1:nrow(Loaloa), 0.7 * nrow(Loaloa))
X_train <- X[train_idx, ]
y_train <- y[train_idx]
m_train <- m[train_idx]
X_test <- X[-train_idx, ]
y_test <- y[-train_idx]
m_test <- m[-train_idx]


p <- y / m # Infection rate
Xbounds <- matrix(c(7.8, 15.3, 3.1, 7.0), ncol = 2, byrow = TRUE)
df <- data.frame(
  lon = X[,1], lat = X[,2], p = p, n = m,
  set = ifelse(1:nrow(Loaloa) %in% train_idx, "Train", "Test")
)

# Obtain African map data (sf object)
africa <- ne_countries(continent = "africa", scale = "medium", returnclass = "sf")

pmap <- ggplot() +
  geom_sf(data = africa, fill = "gray95", color = "gray60") +
  geom_point(data = df,
             aes(x = lon, y = lat, color = p, size = n, shape = set),
             alpha = 0.8) +
  scale_color_viridis_c(name = "Proportion y/m") +
  scale_size_continuous(name = "Trial count m") +
  scale_shape_manual(name = "Dataset", values = c("Train" = 16, "Test" = 17)) +
  coord_sf(xlim = Xbounds[1, ], ylim = Xbounds[2, ], expand = FALSE) +
  theme_minimal() +
  # theme(aspect.ratio = 1) +
  labs(title = "Loaloa infection proportion (y/m): Train vs Test",
       x = "Longitude", y = "Latitude")
ggsave("Loaloa_map.pdf", plot = pmap, width = 10, height = 5)


# Model fitting
Loaloa_bkp_model <- fit_BKP(
  X_train, y_train, m_train, Xbounds, loss = "brier",
  prior = "adaptive", r0 = mean(m_train))

# Print the model summary
summary(Loaloa_bkp_model)

# Plot the fitted model
pdf("Loaloa_bkp.pdf", width = 12, height = 12)
plot(Loaloa_bkp_model)
dev.off()

# LGP model for comparison
Loaloa_gp_model <- gp_init(cf = cf_sexp(), lik = lik_binomial())
Loaloa_gp_model <- gp_optim(
  Loaloa_gp_model, X_train, y_train, trials = m_train,
  method = method_full(), approx = approx_ep(), verbose = FALSE)

# Plot the fitted LGP model
pdf("Loaloa_gp.pdf", width = 12, height = 12)
# Generate the grid for predictions
x1_seq <- seq(Xbounds[1,1], Xbounds[1,2], length.out = 80)
x2_seq <- seq(Xbounds[2,1], Xbounds[2,2], length.out = 80)
grid <- expand.grid(x1 = x1_seq, x2 = x2_seq)

# Predict the grid for plot
prediction <- gp_pred(Loaloa_gp_model, as.matrix(grid),
                      transform = TRUE, var = TRUE, quantiles = c(0.025,0.975))
df <- data.frame(x1 = grid$x1, x2 = grid$x2,
                 Mean = prediction$mean,
                 Upper = prediction$quantiles[,2],
                 Lower = prediction$quantiles[,1],
                 Variance = prediction$var)
dims <- c(1,2)
# Width = prediction$upper - prediction$lower)
p1 <- BKP:::my_2D_plot_fun("Mean", "Predictive Mean", df, dims= dims)
p3 <- BKP:::my_2D_plot_fun("Variance", "Predictive Variance", df, dims= dims)
p2 <- BKP:::my_2D_plot_fun("Upper", "95% CI Upper", df, dims= dims)
p4 <- BKP:::my_2D_plot_fun("Lower", "95% CI Lower", df, dims= dims)
# Arrange into 2×2 layout
grid.arrange(p1, p2, p3, p4, ncol = 2)
dev.off()

# Predict the mean infection rate for the test data
predict_on_test_bkp <- predict(Loaloa_bkp_model, Xnew = X_test)
predict_on_test_gp <- gp_pred(Loaloa_gp_model, xnew = X_test,
                              transform = TRUE, var = TRUE)

# Empirical success rate for the test data
pi_tilde_test <- y_test / m_test

# Mean Squared Error (Brier Score)
mse_bkp <- mean((predict_on_test_bkp$mean - pi_tilde_test)^2)
mse_gp <- mean((predict_on_test_gp$mean - pi_tilde_test)^2)
```

This section presents the workflow of the \pkg{BKP} package based on a real data analysis, including data preprocessing, model fitting, prediction, and visualization. %As an illustrative application, we analyze the prevalence of Loa loa parasite infection in North Cameroon using the \code{Loaloa} dataset \citep{Diggle2007}.  

## Data Preparation 

The \code{Loaloa} dataset, available in the \pkg{mdhglm} package \citep{Lee2018mdhglm}, contains 197 observations collected in North Cameroon between 1991 and 2001.  
For each surveyed location, the dataset records geographical coordinates (\code{longitude}, \code{latitude}), the number of infected individuals (\code{npos}), the total number of surveyed individuals (\code{ntot}), and additional covariates such as altitude and the maximum normalized difference vegetation index (NDVI). 
 
In this analysis, we focus on geographical coordinates as predictors (\code{X}), the number of infected individuals as the response (\code{y}), and the total number surveyed as the trial counts (\code{m}).
The dataset is randomly divided into a training set (70\%) and a testing set (30\%):  
\begin{CodeChunk}
\begin{CodeInput}   
R> data("Loaloa")  
R> # Extract input variables (X), response variable (y), and trial counts (m)
R> X <- as.matrix(Loaloa[, 1:2]) 
R> y <- Loaloa$npos
R> m <- Loaloa$ntot 
R> # Randomly split into training (70%) and testing (30%) sets
R> train_idx <- sample(1:nrow(Loaloa), 0.7 * nrow(Loaloa))
R> X_train <- X[train_idx, ]
R> y_train <- y[train_idx]
R> m_train <- m[train_idx]
R> X_test <- X[-train_idx, ]
R> y_test <- y[-train_idx]
R> m_test <- m[-train_idx] 
\end{CodeInput} 
\end{CodeChunk}   

To investigate the spatial heterogeneity of infection prevalence, Figure \ref{fig:Loaloa_map} presents the observed infection proportion ($y/m$) at each survey location. 
The symbol size is proportional to the sample size, reflecting the relative precision of the observed proportions, while the color gradient encodes the infection rate ranging from 0.0 to 0.5. 
The map further distinguishes between the training set (70\% of locations) and the testing set (30\% of locations), which are used for model fitting and predictive validation, respectively. 
Distinct spatial patterns emerge, with clusters of elevated infection prevalence observed between $10^{\circ}$E and $12^{\circ}$E longitude and another group located between $14^{\circ}$E and $15^{\circ}$E longitude.  


\begin{figure}[!t]
	\centering 
	\includegraphics[width=\linewidth]{Loaloa_map.pdf}  
	\caption{Observed infection proportions ($y/m$) across North Cameroon. Symbol size is proportional to the trial count, color indicates the observed infection rate, and shape distinguishes between training (circle, 70\% of locations) and testing (triangle, 30\% of locations) sets.}
	\label{fig:Loaloa_map}
\end{figure}


## Model fitting and visualization 

We fit the BKP model to the training data using the \fct{fit\_BKP} function.
The model was trained with the Brier score as the loss function, and an adaptive prior was adopted to stabilize estimation in regions with small sample sizes.
The global precision parameter \code{r0} is set to the mean number of trials per location, representing a moderately informative prior that provides smoothing while still allowing the data to contribute substantially to posterior inference.

\begin{CodeChunk}
\begin{CodeInput}
R> Xbounds <- matrix(c(7.8, 15.3, 3.1, 7.0), ncol = 2, byrow = TRUE)
R> Loaloa_bkp_model <- fit_BKP(
+    X_train, y_train, m_train, Xbounds, loss = "brier",
+    prior = "adaptive", r0 = mean(m_train))
R> summary(Loaloa_bkp_model)
\end{CodeInput}
\begin{CodeOutput}

       Beta Kernel Process (BKP) Model   

Number of observations (n):  137
Input dimensionality (d):    2
Kernel type:                 gaussian
Optimized kernel parameters: 0.0076, 0.0460
Minimum achieved loss:       0.00780
Loss function:               brier
Prior type:                  adaptive
    r0: 135.533

Posterior predictive summary (training points):
                      Mean Median     SD Min    Max
Posterior means     0.1508 0.1220 0.1227   0 0.4619
Posterior variances 0.0002 0.0002 0.0003   0 0.0013
\end{CodeOutput}
\end{CodeChunk}  

\begin{CodeChunk} 
\begin{CodeInput}  
R> plot(Loaloa_bkp_model)
\end{CodeInput} 
\end{CodeChunk} 

\begin{figure}[!t]
	\centering 
	\includegraphics[width=\linewidth]{Loaloa_bkp.pdf}  
	\caption{Posterior summaries of the BKP model fitted to the Loa loa dataset, including the predictive mean, variance, and 95\% credible interval surfaces.}
	\label{fig:Loaloa_bkp}
\end{figure}

Figure \ref{fig:Loaloa_bkp} displays the posterior predictive mean, variance, and 95\% credible interval surfaces, providing a comprehensive view of the model's predictive uncertainty.
The posterior mean surface highlights clear spatial variation in Loa loa infection risk, with elevated prevalence predicted in the central–eastern part of the study region. Notably, the predictive mean map from the BKP model effectively captures the major high-prevalence clusters, correctly identifying concentrated yellow and orange regions that correspond to the high-proportion areas observed in Figure \ref{fig:Loaloa_map}. 
The variance surface shows that predictive uncertainty is generally low in areas with dense observations, but increases substantially in peripheral regions where data are sparse.
The 95\% credible interval surfaces further illustrate this pattern: the upper bound suggests that high infection rates cannot be ruled out in some regions, while the lower bound remains close to zero across much of the area.
%Overall, the BKP model not only captures spatial heterogeneity in infection prevalence but also accurately reflects the spatial distribution and concentration of Loa loa infections while quantifying the associated predictive uncertainty.  

For comparison, we also fit a logistic Gaussian process (LGP) model using the \fct{gp\_fit} function from the \pkg{gplite} package \citep{Piironen2022gplite}.  
\begin{CodeChunk}
\begin{CodeInput}   
R> Loaloa_gp_model <- gp_init(cf = cf_sexp(), lik = lik_binomial())
R> Loaloa_gp_model <- gp_optim(
+    Loaloa_gp_model, X_train, y_train, trials = m_train,
+    method = method_full(), approx = approx_ep(), verbose = FALSE)
\end{CodeInput} 
\end{CodeChunk}  

\begin{figure}[!t]
	\centering 
	\includegraphics[width=\linewidth]{Loaloa_gp.pdf}  
	\caption{Posterior summaries of the LGP model fitted to the Loa loa dataset, showing the predictive mean and variance surfaces for comparison with the BKP results.}
	\label{fig:Loaloa_gp}
\end{figure} 

Figure \ref{fig:Loaloa_gp} presents the predictive results of the LGP model. 
Compared with the BKP model, several important differences emerge.  
First, the BKP predictive mean surface is smoother and spatially more coherent, identifying continuous regions of high and low infection prevalence. In contrast, the LGP predictive mean appears more fragmented, with multiple isolated pockets of elevated infection rates. 
Second, the predictive variance differs remarkably: the BKP variance surface is relatively uniform, indicating a consistent level of uncertainty across the study region, whereas the LGP variance surface exhibits highly localized areas of elevated uncertainty (dark purple regions), suggesting less stable predictions in those locations. 
Overall, the BKP model provides a more robust and spatially coherent characterization of Loa loa infection prevalence compared to the LGP model.  

## Predictive performance on test data  
Finally, we evaluate predictive accuracy using predictions from the fitted BKP and LGP models, with performance measured by the Brier score.

\begin{CodeChunk}
\begin{CodeInput}    
R> predict_on_test_bkp <- predict(Loaloa_bkp_model, Xnew = X_test) 
R> predict_on_test_gp <- gp_pred(Loaloa_gp_model, xnew = X_test, 
			transform = TRUE, var = TRUE) 
R> # Calculate the empirical success rate for the test data
R> pi_tilde_test <- y_test / m_test 
R> # Calculate the Mean Squared Error (Brier Score)
R> mse_bkp <- mean((predict_on_test_bkp$mean - pi_tilde_test)^2)
R> mse_gp <- mean((predict_on_test_gp$mean - pi_tilde_test)^2)   
\end{CodeInput}
\end{CodeChunk} 

The results indicate that the BKP model achieves a lower Brier score (0.006) compared to the LGP model (0.019), demonstrating its superior predictive accuracy.  
This shows that the BKP predictions are closer to the observed Loa loa infection rates, and its probabilistic forecasts are therefore more reliable.  

# Summary and discussion {#sec:summary}

This article has presented \pkg{BKP}, a user-friendly \proglang{R} package that implements the Beta Kernel Process (BKP) -- a scalable, interpretable, and fully Bayesian nonparametric framework for modeling spatially varying binomial probabilities. 
The package provides a flexible and modular interface for fitting BKP models to both binary and aggregated binomial data, and extends naturally to the Dirichlet Kernel Process (DKP) for handling multinomial and compositional responses with multiple categories.
To our knowledge, \pkg{BKP} is the first publicly available software for implementing BKP methodology, thereby filling an important gap in the toolkit for spatially varying binomial and multinomial modeling. 

Future development directions include extending the BKP framework to support more complex data structures, such as multivariate responses, functional data, time series, and combinations of qualitative and quantitative covariates. Another promising avenue is to generalize the BKP methodology under alternative likelihoods beyond the binomial family. For example, negative binomial likelihoods are particularly suitable for over-dispersed count data, where the variance exceeds the mean, a common phenomenon in ecological surveys, RNA-seq gene expression counts, and epidemiological incidence data. Geometric likelihoods, as a special case of the negative binomial, naturally model the number of trials until the first success and are useful in reliability analysis, survival studies, and modeling waiting times in event processes. %\url{https://www.datacamp.com/tutorial/negative-binomial-distribution}
These methodological and computational extensions would substantially broaden the applicability of BKP across applied statistics, biostatistics, and machine learning.

At last, we welcome contributions from the community and invite developers to participate in the ongoing maintenance and extension of the package by submitting pull requests via GitHub at \url{https://github.com/Jiangyan-Zhao/BKP/pulls}. 
