how to write the joint density of two correlated uniform random variables? Suppose the marginal distributions of $X$ and $Y$ are both $U[a,b]$, where $U[a,b]$ denotes uniform distribution on $[a,b]$. When $X$ and $Y$ are correlated, how to write their joint density for some prespecified measure of correlation?
 A: I think the right place to start is what correlation means in a nicer setting: Gaussian, of course.
In the bivariate Gaussian distribution, the population correlation is given by the $\rho$ parameter. This tells us something about the fit to a line.
I have tried fitting marginal uniform distributions to a line and failed, but we can still get at the $\rho$ parameter of the bivariate Gaussian.
That $\rho$ parameter is the parameter of the bivariate Gaussian copula, and we can use a Gaussian copula with more than just Gaussian marginal distributions. We can use uniform marginal distributions and specify the parameter in the Gaussian copula, giving us a way to describe the relationship between the two marginal distributions.
When I have simulated this, I have not gotten the sample correlation to be the specified parameter of the Gaussian copula, so I do not know if this quite counts as specifying the correlation. However, the estimate in R's cor function is biased, perhaps accounting for this discrepency. I welcome discussion in the comments about this topic. 
(Remember that a biased estimator isn't inherently a bad estimator; in fact, $s$ is a biased estimator for standard deviation, even though $s^2$ is unbiased for variance, so we use biased estimators all the time.)
library(copula)
set.seed(2020)
N <- 1000
R <- 10000

# define the copula
#
nc <- normalCopula(param = 0.81)

# Define the population distribution with the nc copula and U(0,2) marginals
# 
unif_unif <- mvdc(nc, c("unif","unif"),list(list(min=0, max=2),list(min=0, max=2)))
v <- rep(NA, R)
for (i in 1:R){

  # Sample from the population
  #
  D_uu <- rMvdc(N, unif_unif)

  # Calculate the correlation
  #
  v[i] <- cor(D_uu[,1], D_uu[,2])
  if (i %% 250 == 0){print(i)}
}
plot(density(v))
abline(v=0.81)
mean(v) # 0.7963396, slightly lower than the specified 0.81

EDIT
Let's look at some pictures. Here is the last distribution produced in the simulation: plot(D_uu).

Now compare to a bivariate normal distribution with population correlation $0.81$: plot(qnorm(D_uu[, 1]/2), qnorm(D_uu[, 2]/2)). (This uses the same (Gaussian) copula to relate the marginals but gives each marginal a normal distribution.)

Both plots have the data kind of hugging the line $y=x$, and both see a bulge in the perpendicular direction. However, the behavior of the marginal normal distributions causes more points to be near the mean than in the marginal uniform distributions (consider the PDFs).
