# How to evaluate multivariate log normal distribution at arbitrary point?

**URL:** https://discourse.julialang.org/t/how-to-evaluate-multivariate-log-normal-distribution-at-arbitrary-point/112608
**Category:** General Usage
**Tags:** distributions, monte-carlo, integral, random
**Created:** [April 6, 2024, 2:43pm UTC](https://discourse.julialang.org/t/how-to-evaluate-multivariate-log-normal-distribution-at-arbitrary-point/112608 "2024-04-06T14:43:07Z")
**Posts on this page:** 3
**Page:** 1

<div class="post-metadata">

### Author: ![structural](https://avatars.discourse-cdn.com/v4/letter/s/34f0e0/32.png) [@structural](https://discourse.julialang.org/u/structural)
#### Post date: [April 6, 2024, 2:43pm UTC](https://discourse.julialang.org/t/how-to-evaluate-multivariate-log-normal-distribution-at-arbitrary-point/112608/1 "2024-04-06T14:43:07Z")

</div>

Suppose I have the following distribution over two variables

```julia
dist = MvLogNormal(mean_a_b, cov_a_b)

```

Suppose `h` is an array of dimension 1000\times 2. For each row of array `h`, the first entry represents the rank of the first variable, and the second entry is the rank of the second variable. Is there a way to evaluate the joint distribution `dist` for each of the 1000 points in `h`?

---

<div class="post-metadata">

### Author: ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [April 7, 2024, 5:44am UTC](https://discourse.julialang.org/t/how-to-evaluate-multivariate-log-normal-distribution-at-arbitrary-point/112608/2 "2024-04-07T05:44:11Z")

</div>

The natural thing you can do is use `logpdf` in a `for` loop over each row of `h`. In Julia, this will be essentially as fast as possible. But if for some reason you need vectorization (e.g. to make autodiff faster), you can leverage the fact that Distributions.jl knows a matrix can be a sequence of vector observations. The thing to remember is that each observation needs to be a matrix column, hence the transpose:

```julia
julia> using Distributions, LinearAlgebra

julia> dist = MvLogNormal(zeros(2), I);

julia> h = rand(1000, 2);

julia> logpdf(dist, transpose(h))
1000-element Vector{Float64}:
...

```

If Distributions.jl were less clever, you could also vectorize it manually by wrapping the first argument in a `Ref`, which tells Julia to only broadcast on the second element:

```julia
julia> logpdf.(Ref(dist), eachrow(h))
1000-element Vector{Float64}:
...

```

---

<div class="post-metadata">

### Author: ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [April 7, 2024, 9:46am UTC](https://discourse.julialang.org/t/how-to-evaluate-multivariate-log-normal-distribution-at-arbitrary-point/112608/3 "2024-04-07T09:46:47Z")

</div>

If I understood correctly, the `h` matrix contains _ranks_ of the components in the column of 1000 pairs. This makes the overall distribution dependent across all pairs.  
If it was one column, sampling could be done by sampling, sorting and permuting. If it where several columns and diagonal covariance matrix, the problem again could be split into one column versions. The way it is presented might require some MCMC sampling.
