In [1]:
using Distributions
using PDMats   # this packages is specific to positive definite matrices
$$\forall x \ne 0,\ x^T A x > 0$$
In [2]:
L = [1.0 0 0 ; -2 2 0 ; 1 2 1]
C = L*L'
Out[2]:
3×3 Array{Float64,2}:
  1.0  -2.0  1.0
 -2.0   8.0  2.0
  1.0   2.0  6.0
In [3]:
eig(C)
Out[3]:
([0.0771475,5.50588,9.41697],
[-0.920113 0.36104 0.151799; -0.296781 -0.389819 -0.871758; 0.255565 0.847167 -0.465827])
In [4]:
U = cholfact(C)
Out[4]:
Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:
[1.0 -2.0 1.0; 0.0 2.0 2.0; 0.0 0.0 1.0]
In [5]:
U = [1.0 -2.0 1.0; 0.0 2.0 2.0; 0.0 0.0 1.0]
Out[5]:
3×3 Array{Float64,2}:
 1.0  -2.0  1.0
 0.0   2.0  2.0
 0.0   0.0  1.0
In [6]:
U'*U
Out[6]:
3×3 Array{Float64,2}:
  1.0  -2.0  1.0
 -2.0   8.0  2.0
  1.0   2.0  6.0
In [7]:
cholfact(C,:L)
Out[7]:
Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:
[1.0 0.0 0.0; -2.0 2.0 0.0; 1.0 2.0 1.0]
In [8]:
U'
Out[8]:
3×3 Array{Float64,2}:
  1.0  0.0  0.0
 -2.0  2.0  0.0
  1.0  2.0  1.0
In [9]:
L = [1.0 0.0 0.0; -2.0 2.0 0.0; 1.0 2.0 1.0]
Out[9]:
3×3 Array{Float64,2}:
  1.0  0.0  0.0
 -2.0  2.0  0.0
  1.0  2.0  1.0
In [10]:
μ = [1.0, 2, 1]
Out[10]:
3-element Array{Float64,1}:
 1.0
 2.0
 1.0
In [11]:
Σ = PDMat(C)
Out[11]:
PDMats.PDMat{Float64,Array{Float64,2}}(3,[1.0 -2.0 1.0; -2.0 8.0 2.0; 1.0 2.0 6.0],Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:
[1.0 -2.0 1.0; 0.0 2.0 2.0; 0.0 0.0 1.0])
In [12]:
variance = diag(Σ)
Out[12]:
3-element Array{Float64,1}:
 1.0
 8.0
 6.0
In [13]:
σ = sqrt(variance)
Out[13]:
3-element Array{Float64,1}:
 1.0    
 2.82843
 2.44949
In [14]:
d = MvNormal(μ, C)
Out[14]:
FullNormal(
dim: 3
μ: [1.0,2.0,1.0]
Σ: [1.0 -2.0 1.0; -2.0 8.0 2.0; 1.0 2.0 6.0]
)
In [15]:
d = MvNormal(μ, Σ)
Out[15]:
FullNormal(
dim: 3
μ: [1.0,2.0,1.0]
Σ: [1.0 -2.0 1.0; -2.0 8.0 2.0; 1.0 2.0 6.0]
)
In [16]:
rand(d)
Out[16]:
3-element Array{Float64,1}:
  2.70171
 -3.18036
  3.24845
In [17]:
using RandomStreams

const SEED = 12345

seeds = [SEED, SEED, SEED, SEED, SEED, SEED]
gen = MRG32k3aGen(seeds)
random = next_stream(gen)
Out[17]:
Full state of MRG32k3a generator:
Cg = [12345,12345,12345,12345,12345,12345]
Bg = [12345,12345,12345,12345,12345,12345]
Ig = [12345,12345,12345,12345,12345,12345]
In [18]:
rand_dist(rng::MRG32k3a, Dist::Distribution) = quantile(Dist, rand(rng))
Out[18]:
rand_dist (generic function with 1 method)
In [19]:
d = Normal()
x = Array(Float64, length(μ))

for i =1:length(μ)
    x[i] = rand_dist(random, d)
end
In [20]:
draw = μ + L*x
Out[20]:
3-element Array{Float64,1}:
 -0.140634
  3.33763 
 -1.58243
In [21]:
L = cholfact(C,:L)
Out[21]:
Base.LinAlg.Cholesky{Float64,Array{Float64,2}} with factor:
[1.0 0.0 0.0; -2.0 2.0 0.0; 1.0 2.0 1.0]
In [22]:
draw = μ + L*x

LoadError: MethodError: no method matching *(::Base.LinAlg.Cholesky{Float64,Array{Float64,2}}, ::Array{Float64,1})
Closest candidates are:
  *(::Any, ::Any, !Matched::Any, !Matched::Any...) at operators.jl:138
  *{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},S}(!Matched::Union{Base.ReshapedArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N}}},DenseArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2},SubArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:Union{Base.ReshapedArray{T,N,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N}}},DenseArray},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex,Colon,Int64,Range{Int64}},N}},L}}, ::Union{Base.ReshapedArray{S,1,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N}}},DenseArray{S,1},SubArray{S,1,A<:Union{Base.ReshapedArray{T,N,A<:DenseArray,MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N}}},DenseArray},I<:Tuple{Vararg{Union{Base.AbstractCartesianIndex,Colon,Int64,Range{Int64}},N}},L}}) at linalg/matmul.jl:79
  *(!Matched::Base.LinAlg.AbstractTriangular{T,S<:AbstractArray{T,2}}, ::AbstractArray{T,1}) at linalg/triangular.jl:1496
  ...
while loading In[22], in expression starting on line 1

 in execute_request(::ZMQ.Socket, ::IJulia.Msg) at /home/bastin/.julia/v0.5/IJulia/src/execute_request.jl:169
 in eventloop(::ZMQ.Socket) at /home/bastin/.julia/v0.5/IJulia/src/eventloop.jl:8
 in (::IJulia.##9#15)() at ./task.jl:360
In [ ]: