Bivariate Normal Distribution

n = 9;
(* X is (X1, X2, X3, ..., Xn)' p.49 *)
X = Transpose[Table[RandomVariate[NormalDistribution[1,3],n], {i, 1, 2}]];
TableForm[X,TableHeadings->{{"Trial 1","Trial 2", "Trial 3","Trial 4","Trial 5", "Trial 6", "Trial 7","Trial 8", "Trial 9"},{"X1","X2"}}]

n = 10^3;
σ=3;
μ = 1;
X = Transpose[Table[RandomVariate[NormalDistribution[μ,σ],n], {i, 1, 2}]];

Mean[X]
covx=Covariance[X];
covx//MatrixForm

{Histogram[X[[All,1]]],Histogram[X[[All,2]]]}

Symbolic Computation

σ=.;
(Subscript[σ, 1])=.;
(Subscript[σ, 2])=.;
Σ = .;
ρ  =.;
Σ = {{((Subscript[σ, 1])^2),ρ (Subscript[σ, 1]) (Subscript[σ, 2])},{ρ (Subscript[σ, 1]) (Subscript[σ, 2]),((Subscript[σ, 2])^2)}};
Σ//MatrixForm
PDF[MultinormalDistribution[{(Subscript[μ, 1]),(Subscript[μ, 2])},{{((Subscript[σ, 1])^2),ρ (Subscript[σ, 1]) (Subscript[σ, 2])},{ρ (Subscript[σ, 1]) (Subscript[σ, 2]),((Subscript[σ, 2])^2)}}],{x1,x2}]

Covariance[MultinormalDistribution[{(Subscript[μ, 1]),(Subscript[μ, 2])},Σ]]//MatrixForm

Let’s manually calculate the PDF:

X = .;
n = 2;
x1=.;
x2 = .;
μ=.;
μ1 = .;
μ2 = .;
Σ = {{((Subscript[σ, 1])^2),ρ (Subscript[σ, 1]) (Subscript[σ, 2])},{ρ (Subscript[σ, 1]) (Subscript[σ, 2]),((Subscript[σ, 2])^2)}};
μ = {μ1,μ2};
X = {x1, x2};

pdf  = Simplify[((2*π)^(-n/2))*((Det[Σ])^(-1/2))]*E^Simplify[(-1/2)*((X-μ).Inverse[Σ]).(X-μ)]

(Subscript[σ, 1])=1;
(Subscript[σ, 2]) = 0.5;
ρ = 0.7;
Σ = {{((Subscript[σ, 1])^2),ρ (Subscript[σ, 1]) (Subscript[σ, 2])},{ρ (Subscript[σ, 1]) (Subscript[σ, 2]),((Subscript[σ, 2])^2)}}

data=RandomVariate[d=MultinormalDistribution[{1,2},Σ],10^4];

{Histogram3D[data, ImageSize->300],Plot3D[PDF[d,{x,y}],{x,-4,5},{y,-2,5},PlotPoints->35,PlotRange->All, ImageSize->300]}