distribution = BinormalDistribution[{0,0},{1,(1/2)},(2/3)];

plot = DensityPlot[PDF[distribution][{x , y}] , {x , -3 , 3} , {y , -3 , 3} , ColorFunction->GrayLevel,PlotPoints->50];

(*2D FGP(x , y), UWAGA: pojawiają się kawałki z xy:*)

PDF[distribution][{x,y}]

(*Rysunek, podobnie jak poprzednio:*)

Show[{plot , Graphics[{Blue , PointSize[Tiny],Point/@RandomVariate[distribution , 10000]}]} , FrameLabel->{"X" , "Y"}]

(*normalizacja:*)

Integrate[PDF[distribution][{x , y}] , {x , -∞ , ∞} , {y , -∞ , ∞}]

(*E(X):*)

ex = Integrate[x PDF[distribution][{x , y}] , {x , -∞ , ∞} , {y , -∞ , ∞}]

ey = Integrate[y PDF[distribution][{x , y}] , {x , -∞ , ∞} , {y , -∞ , ∞}]

(*cov(X , Y):*)

covxy = Integrate[(x - ex)(y - ey) PDF[distribution][{x , y}] , {x , -∞ , ∞} , {y , -∞ , ∞}]

(*cov(X , X) = var(X)*)

covxx = Integrate[(x - ex)(x - ex) PDF[distribution][{x , y}] , {x , -∞ , ∞} , {y , -∞ , ∞}]

(*cov(Y , Y) = var(Y)*)

covyy = Integrate[(y - ey)(y - ey) PDF[distribution][{x , y}] , {x , -∞ , ∞} , {y , -∞ , ∞}]

(*macierz kowariancji:*)

{{covxx , covxy} , {covxy , covyy}}//MatrixForm

Covariance[distribution]//MatrixForm

nData = 100000;

data = RandomVariate[distribution , 100000];

exE = Mean[data[[;;,1]]]

eyE = Mean[data[[;;,2]]]

covxyE=Mean[(#[[1]]-exE)(#[[2]]-eyE)&/@data]

covxxE = Mean[(#[[1]]-exE)(#[[1]]-exE)&/@data]

covyyE = Mean[(#[[2]]-eyE)(#[[2]]-eyE)&/@data]

(*Wyliczona "eksperymentalnie" macierz kowariancji (niezbyt porządnym sposobem, lepszy poznamy na następnym wykładzie):*)

{{covxxE , covxyE} , {covxyE , covyyE}}//MatrixForm

Covariance[distribution]//N//MatrixForm

(*czary mary z wektorami własnymi macierzy kowariancji :-)*)

r = (Normalize/@Eigensystem[Covariance[distribution]][[2]]);

(*unitarna macierz transformacji zmiennych losowych:*)

(*{x' , y'} = r.{x , y}*)

Transpose[r]//MatrixForm

r.Transpose[r]//FullSimplify//MatrixForm

(*transformacja macierzy kowariancji jest diagonalna!!!*)

r.Covariance[distribution].Transpose[r]//FullSimplify//MatrixForm

(*rachunek macierzy kowariancji dla nowych zmiennych losowych:*)

data1 = r.#&/@data;

exE1 = Mean[data1[[;;,1]]]

eyE1 = Mean[data1[[;;,2]]]

covxyE1=Mean[(#[[1]]-exE1)(#[[2]]-eyE1)&/@data1]

covxxE1 = Mean[(#[[1]]-exE1)(#[[1]]-exE1)&/@data1]

covyyE1 = Mean[(#[[2]]-eyE1)(#[[2]]-eyE1)&/@data1]

(*wyliczona eksperymentalnie macierz kowariancji dla nowych zmiennych losowych:*)

{{covxxE1 , covxyE1} , {covxyE1 , covyyE1}}//MatrixForm

r.Covariance[distribution].Transpose[r]//FullSimplify//MatrixForm//N

(*macierz r jest Jacobianem!!!*)

D[(r.{x , y})[[1]] , x] == r[[1,1]]

D[(r.{x , y})[[1]] , y] == r[[1,2]]

D[(r.{x , y})[[2]] , x] == r[[2,1]]

D[(r.{x , y})[[2]] , y] == r[[2,2]]

Det[r]//FullSimplify

(*ponieważ wyznacznik jest 1 możemy od razu napisać FGP dla nowych zmiennych losowych*)

(*znikają kawałki mieszane z xy!!!!*)

((PDF[distribution][Transpose[r].{x1 , y1}])/(Abs[Det[r]]))//FullSimplify

Integrate[% , {x1 , -∞ , ∞} , {y1 , -∞ , ∞}]