First, some packages you want to have. Distributions has standard distributions you might want to use. unless you want to code up your own standard distributions. For plotting, you want to add Winston.



In [1]:

    
Pkg.add("Distributions")
Pkg.add("Winston")
Pkg.add("Gaston")









    



INFO: Nothing to be done.
INFO: Nothing to be done.

We're doing importance sampling, so we want to draw from some undrawable distribution "f" using another distrubition "g" that we know how to draw from.

Set your parameters here (e.g.number of samples, mean vector and covariance matrix) here. "eye(ndim)" creates an ndim x ndim identity matrix.



In [43]:

    
using Distributions
using Winston

ndim = 2
nSamples = 10000

mu = zeros(ndim) #array of ndim zeroes
sig = eye(ndim) * 0.8









    



Warning: Possible conflict in library symbol dtrtri_
Warning: Possible conflict in library symbol dgetri_






    Out[43]:





2x2 Array{Float64,2}:
 0.8  0.0
 0.0  0.8

Alright, now let's construct the distribution we AREN'T able to sample from (fDist), and the distribution we ARE able to sample from (gDist). For simplicity, I'm making fDist a slightly different multivariate normal (which we can draw from), but you can drop in whatever function you want.

We'll sample from gDist (samplesFromG) and evaluate fDist and gDist at these values using f(x) and g(x).



In [58]:

    
fDist = MvNormal(mu,2*sig)
gDist = MvNormal(mu,sig)

samplesFromG = rand(gDist, nSamples)

f(x) = pdf(fDist,x)
g(x) = pdf(gDist,x)

sum(f(samplesFromG)./g(samplesFromG))/nSamples

#h = Histogram(f(samplesFromG)/g(samplesFromG))









    



no method Histogram(Array{Float64,2},)
at In[58]:10



In [ ]:

    
# Under construction

function likelihoodFunction(theta::Array{Float64,2}, sigF::Array{Float64,2})
  ndim = size(theta,1)
  ndata = size(theta,2)
  mu = [ mean(theta[i,:]) for i in 1:ndim ]
  likelihood = (2*pi)^(-ndim/2) * det(sigF)^(-0.5) * exp(-*(*(theta-mu,inv(sigF)),transpose(theta-mu) ))
end