Given bivariate data (X, Y), how to determine a cut-off of X that meets some condition of Y?

Question

I have two distributions, $X$ and $Y$ (shown on the horizontal X axis and vertical Y axis, respectively, see image), that represent different ways of scoring some complex system. For a subset of states of a complex system, I have calculated both a score $x$ using the $X$ method, and a score $y$ using the $Y$ method (i.e. the methods are independent), which produces the plot below. Some more points:

• The method for scoring X is quick but can be inaccurate
• The method for scoring Y is slow to compute but extremely robust
• Larger scores are better
• Each distribution has a different scale (the distributions are correlated).

I would like to use the $X$ method as a proxy for the $Y$ method. Is it possible to say that some score $x_{cut}$ in $X$, will map to a score $Y \ge 4$, so for any $x \ge x_{cut}$, we will be certain that the corresponding $y \ge 4$ , with a confidence of say 99.95%? If so, how would I calculate $x_{cut}$?

My naive thinking is to choose the largest $x$ value found for all $y \lt 4$, and use that as my cut-off. But as I've only measured scores for some states in my complex system, I realise that I may not have encountered the largest $x$ that corresponds to $y \lt 4$.

Answer 1

If you are willing to assume the x value is normally distributed you can use Bayes' theorem to calculate the probability of some y score given a value of x.

Bayes' theorem where $X$ is continuous and $Y$ is discrete is given by,

$P(Y=k|X=x) = \frac{f_{X|Y=k}(x)\cdot P(Y=k)}{f_X(x)}$

where $P(Y=k|X=x)$ is the probability of $Y$ is equal to some score $k$, given the x-score is equal to some value $x$. The density function $f_{X|Y=k}$ is the density of $X$ conditional on the y-score. For this density we can use a density, which is normal conditional on $Y$. The numerator is the "global" density of $X$, also we you can assume a normal distribution.

Given some data it is quite easy estimating the different "global" or unconditional probabilities of the y-score ($P(Y=k)$), it's simply the proportions of the scores in the categories 1,2,3,4 and 5.

Also the denominator is quite easy to estimate. The normal distribution is fully characterized by the parameters $\mu$ and $\sigma$, which are the mean and standard deviation, respectively. The empirical counterparts are just the average and the standard deviation of the x-score.

The tricky part is the conditional normal in the numerator. Here you need a set of means, and standard deviations, one for each value of y (1,2,3,4 and 5). So, you need to group the x-values for each corresponding y-value and compute the mean and standard deviation. These values needs to be saved and you can use them to compute the conditional normal for each value of $Y$.

R-code for estimating and computing the conditional probability function

Here is an R example simulating some values and doing the estimation just described. The function conditional_prob computes the probability of some y-score given a x-score. If you want the probability of say Y>=4 given some x-value you can simply sum the probabilities for y=4 and y=5.

# Simulating x and y scores
n <- 100
mean_x <- 35
sd_x <- 2

x <- rnorm(n, mean=mean_x, sd=sd_x)
y <- as.numeric(cut(x, breaks=5))

plot(x,y)

# Estimating parameters for conditional probability
x_grouped_means <- aggregate(x, list(y), mean)
colnames(x_grouped_means) <- c('y', 'mu_hat')

x_grouped_stds <- aggregate(x, list(y), sd)
colnames(x_grouped_stds) <- c('y', 'sigma_hat')

x_global_mean <- mean(x)
x_gloabal_std <- sd(x)

y_relative_freqs <- as.data.frame(table(factor(y))/length(y))
colnames(y_relative_freqs) <- c('y', 'p_hat')


conditional_prob <- function(x, y, grouped_mu_hat_x, grouped_sigma_hat_x, mu_hat_x, sigma_hat_x, p_hat_y) {
  mu_hat_x_given_y <- grouped_mu_hat_x[grouped_mu_hat_x['y']==y, 'mu_hat']
  sigma_hat_x_given_y <- grouped_sigma_hat_x[grouped_sigma_hat_x['y']==y, 'sigma_hat']
  prob_y <- p_hat_y[p_hat_y['y']==y, 'p_hat']
  dnorm(x, mean=mu_hat_x, sd=sigma_hat_x)*prob_y/dnorm(x, mean=mu_hat_x, sd=sigma_hat_x)
}

Compution the probability of Y=>4 given x=35

p4 <- conditional_prob(35, 4, x_grouped_means, x_grouped_stds, x_global_mean, x_gloabal_std, y_relative_freqs)
p5 <- conditional_prob(35, 5, x_grouped_means, x_grouped_stds, x_global_mean, x_gloabal_std, y_relative_freqs)
p4+p5