I’m new to Mathematica and confused about how random variables work.

Say I have a standard normal random variable $ X\sim \mathcal N(0,1)$ and set $ Y := 2X$ . Then $ \operatorname{Cov}(X,Y) = 2$ . Now if I generate $ N$ -iid. draws for $ X$ , the sample covariance of $ X_1,\dots, X_N$ and $ Y_1,\dots, Y_N$ with $ Y_n := 2X_n$ should be close to $ 2$ .

My question is: how can I compute statistics (e.g. sample covariance) in Mathematica? Naively, I tried

`x:=RandomVariate[NormalDistribution[], 5000] y:=2*x Covariance[x,y] `

but the result is clearly incorrect (usually not close to $ 2$ ), because in actuality (I presume) under this model $ Y$ is *independent* of $ X$ .

What is the most convenient way of achieving the correct result?