# Does it make sense to take a weighted average over cross validation results?

I am using P-splines to estimate an unknown function $$f(x)=y$$ that I am fitting to data $$x,y$$. I am using cross validation to estimate the roughness parameter.

It is important to me to not only estimate the function itself, but also its uncertainty at any given point. I am currently partially achieving this through the propagation of the uncertainties in the data, because each $$y$$ value already has estimated 1-standard-deviation uncertainties $$\sigma_y$$. So I do $$N$$ Monte Carlo realizations of the data, re-fit the splines in each realization, and take an average over those results.

However, this only propagates the random uncertainties in the data, and I also want to propagate the systematic uncertainties of the fit. I was thinking I can possibly achieve this by averaging over the cross-validation results. In particular, I was thinking I can take a weighted standard deviation over the fits with different roughness parameters, where the weights are given by the reciprocal mean square error.

Is there statistical justification for this? Or is there already a different well-known way to accomplish what I want?

• Could you not frame it as a Bayesian problem? You would need to think about what your data generation process is though
– Cryo
Commented Jan 17, 2023 at 13:13
• @Cryo I am wondering that too. One "issue" that I'm coming up with is that I can currently solve this problem very easily using least squares (with the known matrix solution to regularized least squares problem), so I can easily estimate the coefficients belonging to my $1000$ p-splines very quickly. I think if I move to a Bayesian paradigm then I will solve for these values using for example MCMC; might be hard. Also, I don't know what to do about the regularization parameter in a Bayesian problem. I guess I put some prior on it and integrate over it though. Do you have some thoughts here? Commented Jan 17, 2023 at 13:58
• @Cryo By the way, I know very well what my data generation process is. The function is $y_j=\int_j K_j(x) f(x) \; \mathrm{d}x$, where $y$ are my data, $K$ are known functions, and $f$ is the (smooth) function I want to estimate. Commented Jan 17, 2023 at 14:10
• 1000 splines may be a stretch for a single MCMC, but do you actually need this? Splines a piecewise polynomials, so I would expect some degree of localization. In this case you should be able to compute spline estimation in one location, then in another etc. Granted, you would still need to connect the edges, so there may be some degree of iteration required.
– Cryo
Commented Jan 17, 2023 at 22:19
• Alternative would be to parametrize the probability distributions, e.g. by normal distributions, and then fit covariance and means using, e.g. maximum likelihood. You would be fitting more coefficients, but it would be some multiple of what you are doing now.
– Cryo
Commented Jan 17, 2023 at 22:24

I gave some replies in the the comments, so here I will try to work with the suggested data generation process, where:

$$y_j = \int dx\,K_j\left(x\right)\cdot f\left(x\right)$$

I am not sure what the $$\int_j$$ means, domain that depends on the variable? I think this can be incorporated into the $$K_j$$, so let the domain stay general.

The aim is to find $$f$$ and the current method of approximation is splines, which can be expressed as:

$$f\left(x\right)=\sum_i\,a_iq^{(i)}\left(x\right) \quad or\quad f\left(x\right)=\sum_{ik}\,a_{ik} q^{(ik)}\left(x\right)$$

Where $$q^{(i)}\left(x\right)$$ could be $$x^i$$, for series expansion, but could also be something like:

$$q^{(ik)}=\begin{cases}x^i,\,x\in \Omega_k\\ 0,\,otherwise\end{cases}$$

which would accommodate splines ($$\Omega_k$$ is some domain where piecewise polynomial is used). You may need to re-use the same coefficients within different spline functions to get continuity right, this can be dealt with appropriately with a similar formalism, it would still be a sum of functions multiplied by some coefficients, which you want to find. Next, let

$$b_j^{ik}=\int\,dx\,K_j\left(x\right) \,q^{(ik)}\left(x\right)$$

These integrals can be computed analytically or numerically, but it should be a finite one-off as I understand, so should be ok.

The problem is then reduced to

$$y_{j}=\sum_{r} b^{(r)}_j a_r$$

Where $$r$$ runs either over $$i$$ or $${ik}$$ combinations. The problem is thus reduced to linear-algebra type exercise. One easy way to tackle it, with uncertainty, is to let $$a_r$$ be multivariate normal (which is what you would have been assuming if you were fitting with least squares anyway), and use conjugate Bayesian approach, which is fully analytic (for multivariate normal likelihood). All you will need to do is to handle the matrix multiplication by $$b^{(r)}_j$$ (it will transform your covariance matrix). So this way you should be able to handle very many params and still get uncertainty handled.

You could also use SVD to decompose the matrix $$b^{(r)}_j$$ and thus reduce it to problem of univariate normal likelihood, and then apply conjgate approach. This will also offer an opportunity to explicitly set some small eigenvalues to zero, which will serve as a regularization

• Thanks for all of this. Yes the $\int_j$ was simply a typo, please ignore the subscript; it should simply be the Fredholm integral equation. Thanks also for the link to the Normal-inverse-Wishart distribution; that appears to be very relevant and useful. I am wondering a few things. I then have to specify the parameters of the prior distribution, right? It's not entirely clear to me how to select those, although I suppose some experimentation is possible. Also, is it the case that this choice of conjugate prior can emulate both a ridge (Tikhonov) and p-spline penalty? (Just a guess.) Commented Jan 18, 2023 at 13:54
• Not sure about penalties, but one can certainly view priors as a form of regularization. Regarding the NIW distribution. That was my first instinct, so I left it, but I recon you would be better off going for SVD first, this would allow you to reduce your set of parameters to be fitted to a set of independent variables, those you could tackle with simple conjugate single-variable priors (so normal-inverse-chi2).
– Cryo
Commented Jan 18, 2023 at 22:16
• How to specify priors? Best situation is when you roughly know where they should be, but for polynomials it may be difficult. Then one can try to arrive at them via prior predictive. Simplify the problem (e.g. do linear polynomials first, reduce number of points), set very primitive priors and see what the predictive looks like. This will give you a sense of what prior is reasonable, then increase the complexity of the polynomials and number of data points. You will not need to get perfectly right, even approx will do
– Cryo
Commented Jan 18, 2023 at 22:21
• Thanks. My impression from this is that the Bayesian approach to this problem is somewhat complicated and involves making several decisions manually. On the other hand, the frequentist approach is to just cross-validate the roughness parameter from the data itself and call it a day, which seems far simpler (basically plug-and-play). But I guess the frequentist approach doesn't give us the same handle on uncertainties, so perhaps that is the tradeoff. Is my impression correct, or is it just my relative inexperience with these methods? Commented Jan 19, 2023 at 10:00
• In many cases, if you opt for Bayesian approach with flat priors, you get frequentist results. You can handle uncertainties in either approach. Frequentist approaches tend to have a lot of read tests that people just take and use (often without understanding all the assumptions and consequences). Bayesian approaches can be more tricky to handle analytically, but these days probabilistic programming frameworks (tensorflow_probability, numpyro) etc, can crack things easily. ...
– Cryo
Commented Jan 19, 2023 at 14:01