# Estimate correlation between data and data-fit model for variance reduction in Monte Carlo estimate

Say that I want to estimate the mean of a function $f$, $\mathbb{E}[f(X)]$, given some input distribution $x\sim P(x)$. I don't know anython about the form of $f$ except that it is smooth and continuous. I can do this with a conventional Monte Carlo estimation:

$\mathbb{E}[f(X)] \approx \frac{1}{N}\sum_{i=1}^N f(x_i)$.

The control variate technique can be used to decrease that variance of the estimate:

$\mathbb{E}[f(X)] \approx \frac{1}{N} \sum_{i=1}^N \left( f(x_i) + \alpha\left(\hat{f}(x_i) - \mathbb{E}[\hat{f}(x)]\right)\right)$.

The varariance of the estiamte is minimized when $\alpha=-\mathrm{CORR}\left(f(x), \hat{f}(x)\right)\sqrt{\frac{\mathrm{VAR}(f(x))}{\mathrm{VAR}(\hat{f}(x))}}$,

where $\mathrm{CORR}$ is the Peason correlation coefficient and $\mathrm{VAR}$ is the variance operator.

It would be nice to use a data fit model as my $\hat{f}$ function. The problem is that it is not clear to me if it is possible to accurately or reliably estimate the Pearson correlation between the truth, $f(x)$, and my data fit model, $\hat{f}(x)$.

To reiterate, the goal is to estimate the Pearson correlation between the truth and the surrogate, $\mathrm{COR}\left(f(x), \hat{f}(x)\right)$, given the probability density $P(x)$, so that I can confidently decrease the variance of the Monte Carlo approximation of $\mathbb{E}[f(x)]$ using the control variate technique with a data-fit model, using a small number of observed $f$, which were observed with random samples of $P(x)$. The problem is that my surrogate will usually have a correlation close to 1 at the observed points, since these are the points used to make the surrogate, and the correlation will be lower at interpolated points. If my surrogate is exceptionally poor, I could increase the variance of my Monte Carlo estimator. So I need to design the correlation estimation procedure to not over-estimate the correlation. Maybe there is some trade-off between the expected variance reduction and the variance of the variance reduction? Ultimately I am looking for a method to best estimate the correlation and decrease the variance of the Monte Carlo estimator.

The only approach I can think of is something like cross-validation. I can leave one or more observations out, construct a data-fit model with the remaining observations, then compute the correlation between the observed points I left out and the values predicted by the data-fit model. But, cross validation does not seem like a reliable solution to me—overestimation of the correlation could lead to an increase in the variance of the control variate Monte Carlo estimator.

Is there a reliable approach, cross validation or otherwise, for estimating the correlation between a data-fit model and observed data in this setting?

Everything below this line is to illustrate my question. Feel free to skip this part.

Illustrative Example

Here, I implemented a leave-one-out cross-validation approach to estimating the correlation, assuming that the input $x$ is one-dimensional and uniformly distributed from 0 to 10 and $f$ has a smooth and hard to predict form. I used a one-dimensional input to illustrate my point here, but my real problem has a much larger input dimension.

I used a Gaussian Process as the surrogate, though the important aspect of that is that it interpolates across the observed data. For each point in my observations, I left one out and fit the Gaussian Process to the remaining data points, then recorded the surrogate's prediction and the truth value. I decided not to leave the end points out so that I was only using the Gaussian Process for interpolation. I used these pairs of truth and predicted values to estimate the correlation between the function and the surrogate.

This was just my take on the problem. Ultimately, I'd like to experiment with other cross-validating approaches. I'm not sure what other approaches I could try besides cross-validation. I'm not sure if cross-validation is even a valid approach to this problem.

Coded implementation

import numpy as np
from matplotlib import pyplot as plt

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C

np.random.seed(5)

# this is the truth. In my real problem, there is no known functional form and I have limited observations.
def f(x):
return x * np.sin(x) - np.cos(5 * x) * np.sqrt(x)

# ----------------------------------------------------------------------
#  assume 20 observed inputs that are uniformly distributed
XL = 0
XU = 10
x_observed = np.atleast_2d(sorted(np.random.uniform(XL, XU, 20)))

# record of predictions
y_pred = []

# Instanciate kernel
kernel = C(1.0, (1e-3, 1e3)) * RBF(10, (1e-2, 1e2))

# estimate correlation between surrogate and truth
#   make the surrogate but leave one out at a time. Always use the edge points.
for ii in range(1, x_observed.shape-1):

# prepare truth data leaving each one out at a time but never excluding the sides
X = x_observed[:, list(range(0, ii)) + list(range(ii+1, x_observed.shape))].T
y = f(X.T).ravel()

# Fit to data using Maximum Likelihood Estimation of the parameters
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)
gp.fit(X, y[:, None])

# Make the prediction over all observations
y_pred.append(gp.predict(x_observed[:, ii:ii+1], return_std=False))

# Create surrogate using all observed data points
X = x_observed.T
y = f(X.T).ravel()
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)
gp.fit(X, y[:, None])

# create points to estimate true correlation (we are using a uniform distribution)
xs = np.random.uniform(XL, XU, 100000)

# create domain for plotting
xb = np.linspace(XL, XU, 10000)

# estimate correlation
est_corr = np.corrcoef(f(x_observed)[0, 1:-1], y_pred)

# calculate tru correlation
true_corr = np.corrcoef(f(xs), gp.predict(np.atleast_2d(xs).T)[:, 0])

# estimate mean with monte carlo
hf = np.mean(f(x_observed)[0, 1:-1])

# estimate mean with GP using small number of samples
gphf = np.mean(gp.predict(np.atleast_2d(x_observed).T)[:, 0])

# estimate mean with GP using large number of samples
gpmean = np.mean(gp.predict(np.atleast_2d(xs).T)[:, 0])

# estimate the control variate value
alpha_est = -1 * np.corrcoef(f(x_observed)[0, 1:-1], y_pred) * np.std(f(x_observed)[0, 1:-1]) / np.std(gp.predict(np.atleast_2d(x_observed).T)[:, 0])

# calculate the true control variate value
alpha = -1 * np.corrcoef(f(xs), gp.predict(np.atleast_2d(xs).T)[:, -1]) * np.std(f(xs)) / np.std(gp.predict(np.atleast_2d(xs).T)[:, 0])

print("I estimate the correlation of your surrogate is %04.3f" % est_corr)
print("True correlation: %04.3f" % true_corr)
print("True mean %f"%np.mean(f(xs)))
print("Monte Carlo Estimated mean %f"%hf)
print("Conntrol Variate Mean with estimated Correlation %f"% ( hf + alpha_est * (gphf - gpmean) ))
print("Conntrol Variate Mean with actual Correlation %f"% ( hf + alpha * (gphf - gpmean) ))
print("-------")
print("Approach | Error")
print("CV (estimated corr) | %f"%np.abs(np.mean(f(xs)) - (hf + alpha_est * (gphf - gpmean) )))
print("CV (actual corr) | %f"%np.abs(np.mean(f(xs)) - (hf + alpha * (gphf - gpmean) )))
print("Monte Carlo | %f"%np.abs(np.mean(f(xs)) - (hf)))

# Plot the truth, the surrogate, and the observed points
plt.plot(xb, f(xb), label='$f(x)$')
plt.plot(x_observed[0, :], f(x_observed[0, :]), 'r.', markersize=10, label=u'Observed $f(x)$')
plt.plot(xb, gp.predict(np.atleast_2d(xb).T), label=r'$\widehat{f}(x)$')
plt.xlabel('$x$')
plt.legend(loc='upper left')
plt.savefig('fvsfhat')
plt.clf()

# plot the points used to estimate the correlation and the points used to calculate the true correlation
plt.scatter(f(xs), gp.predict(np.atleast_2d(xs).T), c='k', marker='o', s=1, label='Data used to calculate actual corerelation.')
plt.scatter((x_observed)[0, 1:-1], y_pred, c='purple', marker='*', label='Data used to calculate the correlation estimate.\nThese were obtained using leave-one-out "cross-validation."')
plt.xlabel(r'$f(x)$')
plt.ylabel(r'$\widehat{f}(x)$')
plt.legend()
plt.savefig('correlationData')
plt.clf()


Output:

I estimate the correlation of your surrogate is 0.962

True correlation: 0.856

True mean 0.823682

Monte Carlo Estimated mean 0.226071

Conntrol Variate Mean with estimated Correlation 0.698964

Conntrol Variate Mean with actual Correlation 0.720262

Approach | Error

CV (estimated corr) | 0.124718

CV (actual corr) | 0.103420

Monte Carlo | 0.597611

(In the above output, both control variate estimates were more accurate than the conventional Monte Carlo estimator. This trend should hold for most random number generator seeds.)

I want to accurately estimate the correlation between my truth and surrogate given the prescribed input distribution. In the above example, I overestimated the correlation (my estimate, 0.96, was much greater than the actual correlation, which I calculated as 0.86) by leaving one observation out at a time, recording the surrogate's prediction of the out point, then computing the correlation between the predicted points and the truth values.

I compared the accuracy of the control variate Monte Carlo method using the leave-one-out cross-validation estimate of the correlation to a scenario where I somehow knew a very accurate estimate of the correlation, even only with a small sample size available. Here, I used 6 samples to estimate the mean value. The plot compares the absolute errors associated with the control variate Monte Carlo approaches assuming that I have an accurate estimation of the correlation against if I use the leave-one-out cross-validation estimate of the correlation.

import numpy as np
from matplotlib import pyplot as plt

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C

np.random.seed(5)

# this is the truth. In my real problem, there is no known functional form and I have limited observations.
def f(x):
return x * np.sin(x) - np.cos(5 * x) * np.sqrt(x)

# ----------------------------------------------------------------------
#  assume observed inputs that are uniformly distributed
XL = 0
XU = 10
NSAMPS = 6

def experiment():
# make observations
x_observed = np.atleast_2d(sorted(np.random.uniform(XL, XU, NSAMPS)))

# record of predictions
y_pred = []

# Instanciate kernel
kernel = C(1.0, (1e-3, 1e3)) * RBF(10, (1e-2, 1e2))

# estimate correlation between surrogate and truth
#   make the surrogate but leave one out at a time. Always use the edge points.
for ii in range(1, x_observed.shape-1):

# prepare truth data leaving each one out at a time but never excluding the sides
X = x_observed[:, list(range(0, ii)) + list(range(ii+1, x_observed.shape))].T
y = f(X.T).ravel()

# Fit to data using Maximum Likelihood Estimation of the parameters
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)
gp.fit(X, y[:, None])

# Make the prediction over all observations
y_pred.append(gp.predict(x_observed[:, ii:ii+1], return_std=False))

# Create surrogate using all observed data points
X = x_observed.T
y = f(X.T).ravel()
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)
gp.fit(X, y[:, None])

# create points to estimate true correlation (we are using a uniform distribution)
xs = np.random.uniform(XL, XU, 100000)

# create domain for plotting
xb = np.linspace(XL, XU, 10000)

# estimate correlation
est_corr = np.corrcoef(f(x_observed)[0, 1:-1], y_pred)

# calculate tru correlation
true_corr = np.corrcoef(f(xs), gp.predict(np.atleast_2d(xs).T)[:, 0])

# estimate mean with monte carlo
hf = np.mean(f(x_observed)[0, 1:-1])

# estimate mean with GP using small number of samples
gphf = np.mean(gp.predict(np.atleast_2d(x_observed).T)[:, 0])

# estimate mean with GP using large number of samples
gpmean = np.mean(gp.predict(np.atleast_2d(xs).T)[:, 0])

# estimate the control variate value
alpha_est = -1 * np.corrcoef(f(x_observed)[0, 1:-1], y_pred) * np.std(f(x_observed)[0, 1:-1]) / np.std(gp.predict(np.atleast_2d(x_observed).T)[:, 0])

# calculate the true control variate value
alpha = -1 * np.corrcoef(f(xs), gp.predict(np.atleast_2d(xs).T)[:, -1]) * np.std(f(xs)) / np.std(gp.predict(np.atleast_2d(xs).T)[:, 0])

return( ( hf + alpha_est * (gphf - gpmean) ), ( hf + alpha * (gphf - gpmean) ))

exps = [experiment() for _ in range(30)]
x, y = zip(*exps)
plt.scatter(np.abs(x), np.abs(y))
plt.title("Estimated Monte Carlo Estimate Errors.\nThe Estimated means used %i samples."%NSAMPS)
plt.xlabel('Absolute Error, Leave-One-Out Approach')
plt.ylabel('Absolute Error, Using Correlation Esimated With Very Large N')
plt.draw()
plt.savefig('errorPlot')
plt.clf() • Without any a priori knowledge/assumptions the function $\hat{f}(x)$ can not provide any more information about your data sample $f(x)$ and can not reduce the variance of the estimate of the mean. So, if you would have such information, about $f(x)$ then what about using it more directly instead of via the variance reduction technique? – Sextus Empiricus Jul 4 '18 at 10:06

### You are not exactly using the control variate method

What you are basically doing is using a weighted sum with the mean of the surrogate function (*) note that this is not the expectation value, see more in the comment below to estimate the true the mean.

$$\begin{array}\\ \frac{1}{N} \sum_{i=1}^N \left( f(x_i) + \alpha\left(\hat{f}(x_i) - \mathbb{E}[\hat{f}(x)]\right)\right) & = & (1+\alpha) \overline{f(x_i)} - \alpha \mathbb{E}[\hat{f}(x)] \\ &= & (1+\alpha) \overline{f(x_i)} - \alpha \overline{\hat{f}(x)} \end{array}$$

since $f(x_i)=\hat{f}(x_i)$. (Note that your hfand gphf variables are equal, except that hf is a mean for 18 points, without the boundary, and gphf for 20 points)

(*)This is a different situation than the control variate technique for which $\hat{f}(x)$ is a fixed function instead of depending on the sampled $x_i$ and $f(x_i)$. Your computation of $\mathbb{E}(\hat{f}(x))$

gpmean = np.mean(gp.predict(np.atleast_2d(xs).T)[:, 0])

is not a true expectation value. It is the mean of the observation dependent fitted function. This fit can vary and based on that variation you should compute the expectation value. This is why I used $\overline{\hat{f}(x)}$ instead of $\mathbb{E}[\hat{f}(x)]$.

So this is not an unbiased estimator since $\mathbb{E}[\overline{\hat{f}(x)}]$ may possibly be different from $\mathbb{E}[f(x)]$

• answer 1: Your method would equal the control variate method if, after obtaining $\hat{f}(x)$ from a set of $(x_i,f(x_i))$, you draw a new set of $(x_i^\prime,f(x_i^\prime))$, and perform the Monte Carlo estimate with these new data. In this case you can estimate the correlation by comparing $f(x_i^\prime)$ and $\hat{f}(x_i^\prime)$.

### What you should do with your alternative method

Now if the point is to find a value for $\alpha$ which minimizes the variation of the estimator (note this may not be the ultimate goal since it is about the total error which is variance plus bias), then the question is not about the correlation between $f(x_i)$ and $\hat{f}(x_i)$ but about the correlation between the means $m=\overline{f(x_i)}$ and $m^\star = \overline{\hat{f}(x)}$.

The variable

$$(1+\alpha) m - \alpha m^\star$$

will have variance as:

$$(1+\alpha)^2 \text{Var}(m) + \alpha^2 \text{Var}(m^\star) - 2 (1+\alpha)\alpha \text{Cov}(m,m^\star)$$

• answer 2 If you have sufficient data then you might be able to use some cross validation scheme to learn an experimental rule of thumb about typical values for $\rho$ and $\frac{\text{Var}(m^\star)}{\text{Var}(m)}$ and what value to be used for the parameter $\alpha$ (I think it would not be good to use the correlation based on a single sample and function).

If you have more knowledge about the type of functions that you may encounter then you could possibly obtain a theoretically derivation for $\rho$ and $\frac{\text{Var}(m^\star)}{\text{Var}(m)}$ or how to efficiently estimate it experimentally.

• Thanks for your really nice answer. I'm really interested in your last point about having prior knowledge about $f$. Say $f(x)$ is a highly nonlinear partial differential equation. If I had SVD information about $f(x)$, or if I had recorded $\nabla f(x_i)$, could I use that information to somehow inform a procedure to estimate the correlation between $f(x)$ and the data-fit model $\hat{f}$? – kilojoules Jul 4 '18 at 15:51
• Also I don't understand your point about this procedure being biased—in the limit of infinite sample points it seems like the first equation you wrote converge to the mean. – kilojoules Jul 4 '18 at 15:53
• That makes it asymptotically unbiased if the function $\hat{f}$ converges to $f$. But you may have a situation that the bias for a small finite sample is already small or zero. I just mentioned that $E[\overline {\hat {f}}]$ is not necessarily equal to $E[f(x)]$ whereas the expression for the control variate method ensures an unbiased estimator. – Sextus Empiricus Jul 4 '18 at 16:13
• I did a simulation with your function repeating the computation 200 times you get a correlation of ~0.4 and a ratio $var(m^\star)/var(m)$ ~0.3 . Then an $\alpha$ of approximately 0.9 is best. You could do several other tests to make sure this is general in other cases but for reasonable large variations you still get that an $\alpha$ about -0.9 is much better than an $\alpha$ at 0. If in doubt you can alway use an $\alpha$ at -1, which means that you purely use the value of the surrogate and make no use of the possible correlation. – Sextus Empiricus Jul 4 '18 at 16:27
• an open question is whether the control variate method (but using new data after establishing $\hat{f}$) is better than simply using the mean of $\hat{f}$ (or possibly you could use a combination of the two). – Sextus Empiricus Jul 4 '18 at 16:30