  In many applications that involve the inference of an unknown smooth
function, the inference of its derivatives will often be just as important as
that of the function itself. To make joint inferences of the function and its
derivatives, a class of Gaussian processes called $p^{\text{th}}$ order
Integrated Wiener's Process (IWP), is considered. Methods for constructing a
finite element (FEM) approximation of an IWP exist but have focused only on the
order $p = 2$ case which does not allow appropriate inference for derivatives,
and their computational feasibility relies on additional approximation to the
FEM itself. In this article, we propose an alternative FEM approximation,
called overlapping splines (O-spline), which pursues computational feasibility
directly through the choice of test functions, and mirrors the construction of
an IWP as the Ospline results from the multiple integrations of these same test
functions. The O-spline approximation applies for any order $p \in
\mathbb{Z}^+$, is computationally efficient and provides consistent inference
for all derivatives up to order $p-1$. It is shown both theoretically, and
empirically through simulation, that the O-spline approximation converges to
the true IWP as the number of knots increases. We further provide a unified and
interpretable way to define priors for the smoothing parameter based on the
notion of predictive standard deviation (PSD), which is invariant to the order
$p$ and the placement of the knot. Finally, we demonstrate the practical use of
the O-spline approximation through simulation studies and an analysis of COVID
death rates where the inference is carried on both the function and its
derivatives where the latter has an important interpretation in terms of the
course of the pandemic.
