In [23]:
using Distributions, StatsBase, StatsPlots, LaTeXStrings
# Distribution: package with all common distribution
# StatsBase: some common statistics operations
# StatsPlots: Plot + some special stat plots
# LaTeXStrings: Latex writing in figures

Ex 3

In [8]:
n = 1000
X = rand(Cauchy(),n)
# In Julia element wise operations like take the inverse of each element of an array is denoted with a dot . 
# The purpose is to clarify what operation on what object is really done.
Y = 1 ./rand(Cauchy(),n) 
histogram(X,alpha=0.5,label=L"X",xlims=(0,100),nbins=1000)
histogram!(Y,alpha=0.5,label= L"Y",xlims=(0,100),nbins=1000)
Out[8]:

Ex 4

In [6]:
n = 10_000
d = Normal()
X = rand(d,n)
Y = -X
histogram([X,Y],alpha=0.5,label=[L"X" L"Y"])
Out[6]:

Ex 6

In [9]:
BienayméChebychev(d,a) = var(d)/a^2
Cantelli(d,a) = var(d)/(a^2+var(d))
Out[9]:
Cantelli (generic function with 1 method)
In [10]:
d = Normal(0,1)
xr = 0:0.01:5
plot(x-> 1-cdf(d,x),xr, label=L"P(X\geq a)")
plot!(x-> BienayméChebychev(d,x), xr, label="Bienaymé Chebychev")
plot!(x-> Cantelli(d,x), xr, label="Cantelli")
title!("$d")
ylims!(0,1)
xlabel!(L"a")
Out[10]:

Are theses bounds sharp? i.e. is the bound reached by some distributions?

Exemple: $\mathbb{P}(X=-1) = \mathbb{P}(X=-1) = 1/(2k^2)$ and $ \mathbb{P}(X=0) = 1-1/k^2$ for $k\geq 1$

In [49]:
k = 2 # k≥1
xs = [-1.,0,1] # Values of X
ps = [1/(2*k^2),1-1/k^2, 1/(2*k^2)] # Associated probability
d = DiscreteNonParametric(xs, ps)

println(L"$\mathbb{E}(X) = $", mean(d)) # 0
println(L"\operatorname{Var}(X) = ", var(d)) # 1/k^2

xr = 0:0.01:5
plot(x-> 2(1-cdf(d,x)),xr, label=L"P(|X|\geq a)")
plot!(x-> BienayméChebychev(d,x), xr, label="Bienaymé Chebychev")
# plot!(x-> Cantelli(d,x), xr, label="Cantelli")
ylims!(0,1)
xlabel!(L"a")

$\mathbb{E}(X) = $0.0
$\operatorname{Var}(X) = $0.25
Out[49]:

Exemple: $\mathbb{P}(X=1) = \dfrac{v}{1+v}$ and $\mathbb{P}(X=-v) = \dfrac{1}{1+v}$ for $v\geq 0$

In [13]:
v = 0.2;
xs = [1.,-v] # Values of X
ps = [v/(1+v), 1/(1+v)] # Associated probability
d = DiscreteNonParametric(xs, ps)

println(L"\mathbb{E}(X) = ", mean(d))
println(L"\operatorname{Var}(X) = ", var(d))

xr = 0:0.01:5
plot(x-> 1-cdf(d,x),xr, label=L"P(X\geq a)")
# plot!(x-> BienayméChebychev(d,x), xr, label="Bienaymé Chebychev")
plot!(x-> Cantelli(d,x), xr, label="Cantelli",c=3)
ylims!(0,1)
xlabel!("a")

$\mathbb{E}(X) = $0.0
$\operatorname{Var}(X) = $0.2
Out[13]:

Ex 7

In [15]:
f7(x,y) = 2/Ï€*exp(-x*(1+y^2))
xr = range(0, stop = 2, length = 100)
yr = range(0, stop = 3, length = 100)

contourf(xr, yr, (xr, yr) -> f7(xr,yr))
Out[15]:

Ex 8

To simulate $f_{X,Y}(x,y)=(1+xy)\mathbb{1}_{(x,y)\in [-1,1]^2}/4$ one can use the conditional density $f_{X|Y}(x,y) = (1+xy)\mathbb{1}_{x\in [-1,1]}/2$ and marginal $f_{Y}(y) = \mathbb{1}_{y\in [-1,1]}/2$ (= uniform density) so that $$f_{X,Y}(x,y)=f_{X|Y}(x,y)\times f_{Y}(y)$$ In this latter form to simulate a pair $(X,Y)$ one just need to draw $Y\sim U([-1,1])$ and then draw from the conditional density $f_{X|Y}(x,y)$.

We have $$ F_{X|Y}(t) = \mathbb{P}(X<t|Y=y) = \int_{-1}^x \frac{1}{2} (1+x y) \, dx = \frac{1}{4} (1+x) ((x-1) y+2)$$

To use the Inverse Transform method we need to compute $F^{-1}_{X|Y}(U)$ which here can be found as $$F^{-1}_{X|Y}(U) = \frac{\sqrt{4 U y+y^2-2 y+1}-1}{y}$$

Plot the distributions

In [16]:
f(x,y) = (1+x*y)/4
g(x,y) = 1/4
Out[16]:
g (generic function with 1 method)
In [17]:
xr = yr = range(-1, stop = 1, length = 100)
contour(xr, yr, (xr, yr) -> f(xr,yr), fill = true)
Out[17]:
In [18]:
xr = yr = range(-1, stop = 1, length = 100)
contour(xr, yr, (xr, yr) -> g(xr,yr), fill = true)
Out[18]:

Generate samples from $f$ and $g$

In [19]:
# F^{-1}_{X|Y=y}(U)
sampler_x(u,y) = (-1 + sqrt(1 - 2*y + 4*u*y + y^2))/y
Out[19]:
sampler_x (generic function with 1 method)
In [20]:
n = 100_000

x_ind = rand(Uniform(-1,1),n)
y_ind = rand(Uniform(-1,1),n)

y = rand(Uniform(-1,1),n)
x = sampler_x.(rand(n),y)
Out[20]:
100000-element Vector{Float64}:
 -0.04675771816745691
  0.5164990185572046
  0.09743529134896399
  0.21543762646852546
 -0.4816104749486907
  0.053689368138664983
 -0.8974934449742633
 -0.544247942952904
 -0.7478569686154085
  0.1415653600005768
  0.7093047015974552
 -0.4661396629229515
 -0.2101380155805106
  â‹®
  0.10053237380239684
 -0.8049179716442108
  0.3263020704742913
  0.9808425664737404
 -0.9960090667936562
  0.4307726208481279
 -0.37999659447977374
 -0.16164545651688325
  0.2900075269349879
 -0.3659787462653613
  0.01608689849591435
  0.11481559096421551

Plot the simulated densities

Density of $f_{X,Y}$ with level lines (contour)

In [21]:
h = fit(Histogram, (x,y), nbins=10)
contour(midpoints(h.edges[1]),midpoints(h.edges[2]),h.weights,xlabel=L"x", ylabel=L"y", fill = true)
Out[21]:

2D histogram of $f_{X,Y}$ and $g_{X',Y'}$ with the marginal (1D) histogram of $f_{X}$, $f_{Y}$, $g_{X'}$ and $g_{Y'}$

In [22]:
# , title = L"f(x,y)=(1+xy)\mathbb{1}_{(x,y)\in [-1,1]^2}/4"
p1 = marginalhist(x,y,xlabel=L"x", ylabel=L"y",bins=20) 
#, title = L"f(x,y)=\mathbb{1}_{(x,y)\in [-1,1]^2}/4"
p2 = marginalhist(x_ind,y_ind,xlabel=L"x", ylabel=L"y",bins=100)
plot(p1, p2, size=(1000,500))
Out[22]: