using Distributions, PyPlot N = 100 generative_dist = MvNormal([0,1.], [0.8 0.5; 0.5 1.0]) function plotObservations(obs::Matrix) plot(obs[1,:], obs[2,:], "kx", zorder=3) fill_between([0., 2.], 1., 2., color="k", alpha=0.4, zorder=2) # Shaded area text(2.05, 1.8, "S", fontsize=12) xlim([-3,3]); ylim([-2,4]); xlabel("a"); ylabel("b") end D = rand(generative_dist, N) # Generate observations from generative_dist plotObservations(D) x_dot = rand(generative_dist) # Generate x∙ plot(x_dot[1], x_dot[2], "ro"); using HCubature, LinearAlgebra# Numerical integration package # Maximum likelihood estimation of 2D Gaussian μ = 1/N * sum(D,dims=2)[:,1] D_min_μ = D - repeat(μ, 1, N) Σ = Hermitian(1/N * D_min_μ*D_min_μ') m = MvNormal(μ, convert(Matrix, Σ)); # Contour plot of estimated Gaussian density A = Matrix{Float64}(undef,100,100); B = Matrix{Float64}(undef,100,100) density = Matrix{Float64}(undef,100,100) for i=1:100 for j=1:100 A[i,j] = a = (i-1)*6/100 .- 2 B[i,j] = b = (j-1)*6/100 .- 3 density[i,j] = pdf(m, [a,b]) end end c = contour(A, B, density, 6, zorder=1) PyPlot.set_cmap("cool") clabel(c, inline=1, fontsize=10) # Plot observations, x∙, and the countours of the estimated Gausian density plotObservations(D) plot(x_dot[1], x_dot[2], "ro") # Numerical integration of p(x|m) over S: (val,err) = hcubature((x)->pdf(m,x), [0., 1.], [2., 2.]) println("p(x⋅∈S|m) ≈ $(val)") open("../../styles/aipstyle.html") do f display("text/html", read(f, String)) end