Generating data that satisfies certain mean, median and covariance in python I would like to generate pseudo-random data of 16,000 instances with two positive features each. One with a mean value of 6300 and median of 4600 and the other feature with mean of 12500 and median of 12000.
Up until now I managed to get close to what I desire by generating two random batches of below and above the median and then tweaking it a bit to get the required means. The problem now is that the covariance between the two features needs to be equal to 0.97*σ(feature1)*σ(feature2), and I am lost in how to generate the whole data with these requirements. Is there some package or function for generating data with specific values?
 A: For simplicity of the equations, we simplify the problem to deal with means equal to zero.
Then the problem becomes

*

*X with a mean of 0 and median of -1700

*Y with a mean of 0 and median of -500.

A simple way is to sample from a distribution with three points $(X_1,Y_1)$ ,$(X_2,Y_2)$ and $(X_3,Y_3)$ with weights $w_1,w_2,w_3$.
You can find some weights that solve your problem by fiddling around. One solution is
$$w_1 = 0.25, w_2 = 0.5, w_3 = 0.25$$
$$\begin{array}{}
(X_1,Y_1) &=& (400,-147.2911)\\ 
(X_2,Y_2) &=& (-1700,-500)\\
(X_3,Y_3) &=& (3000,1147.291)
\end{array}$$
The strategy is to fix some numbers and then solve the remaining ones by a computer.
One usefull trick is to  set the weight of $w_2 = 0.5$ and $X_2 = -1700$ and $Y_2 = -500$ such that the median is fixed and the remaining equations are easy to solve.
In R the code would be
### function that needs to be 0
f <- function(Y_3) {
  X_3 = 3*10^3
  X_2 = -1700
  Y_2 = -500
  X_1 = -(X_3+2*X_2)
  Y_1 = -(Y_3+2*Y_2)
  X = c(X_1,X_2,X_2,X_3)
  Y = c(Y_1,Y_2,Y_2,Y_3)
  cor(X,Y)-0.97
}
f <- Vectorize(f)


### plot to visually see the maximum and minimum
Y_3 = seq(-10^4,10^4,1)
plot(Y_3, f(Y_3))

### root finding function
Y_3 = uniroot(f, c(832,10^4))$root



### these are the values
X_3 = 3*10^3
Y_3 ### this one is computed by the root finding
X_2 = -1700
Y_2 = -500
X_1 = -(X_3+2*X_2)
Y_1 = -(Y_3+2*Y_2)

The sampling is finally done by sampling the categorical distribution. With weights $w_1 = w_3 = 0.25$ and $w_2 = 0.5$
