# Generating high-dimensional Sobol sequences

**URL:** <https://discourse.julialang.org/t/generating-high-dimensional-sobol-sequences/21521>\
**Category:** Modelling & Simulations\
**Created:** [March 6, 2019, 5:39am UTC](https://discourse.julialang.org/t/generating-high-dimensional-sobol-sequences/21521 "2019-03-06T05:39:43Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![jocklawrie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jocklawrie/32/33616_2.png) [@jocklawrie](https://discourse.julialang.org/u/jocklawrie)\
**Post date:** [March 6, 2019, 5:39am UTC](https://discourse.julialang.org/t/generating-high-dimensional-sobol-sequences/21521/1 "2019-03-06T05:39:43Z")

</div>

Hi there,

I am using Sobol sequences to generate say 10,000 draws from the Multivariate Normal with dimension 100, mean vector 0 and covariance matrix equal to the identity matrix.

Generating the sequence is fine, but I find that the distribution of the resulting points is highly concentrated. That is, the density of each draw, relative to that of the other draws, is almost zero for 99% of points, then several orders of magnitude higher for just a few points. I’d expect the relative densities to be a little more even.

Example below. Any ideas what I’m doing wrong?

Thanks

```julia
using Sobol
using Distributions
using UnicodePlots

function generate_z(K::Int, nr::Int)
    # Initialize result with draws from the Uniform distribution over the unit cube [0, 1]^K.
    result = fill(0.0, K, nr)
    s = SobolSeq(K)
    skip(s, nr)
    for r = 1:nr
        result[:, r] .= next!(s)
    end

    # Convert to draws from MVN(z | 0_k, I_k) using the inverse CDF
    d = Normal(0.0, 1.0)
    for r = 1:nr
        for i = 1:K
            result[i, r] = quantile(d, result[i, r])
        end
    end
    result
end

function zdensity(z::Array{Float64, 2})
    K, nr = size(z)
    result = fill(1.0, nr)
    d = Normal(0.0, 1.0)
    totaldens = 0.0
    for r = 1:nr
        for i = 1:K
            result[r] *= pdf(d, z[i, r])
        end
        totaldens += result[r]
    end
    result ./= totaldens
end

z = generate_z(100, 10_000);
dens = zdensity(z)
lineplot(dens)
describe(dens)

```

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [March 6, 2019, 8:47am UTC](https://discourse.julialang.org/t/generating-high-dimensional-sobol-sequences/21521/2 "2019-03-06T08:47:03Z")

</div>

First, I would use the sum of **log** pdfs, which is much less prone to numerical error. These look fine to me when compared with random normals. Also simple statistics look ok.

```julia
using Statistics, LinearAlgebra
extrema(mapslices(mean, z; dims = 2))
extrema(mapslices(std, z; dims = 2))
z_log_density(z) = normalize!(vec(sum(logpdf.(Ref(Normal(0, 1)), z); dims = 1)), 1)
z = generate_z(100, 10_000)
z_dens = z_log_density(z)
sim_z = randn(size(z))
sim_z_dens = z_log_density(sim_z)
ps = 0:0.1:1
quantile(z_dens, ps)
quantile(sim_z_dens, ps)

```

yields

```julia
julia> extrema(mapslices(mean, z; dims = 2))
(-0.0007172362463843937, 0.0009932045561835572)

julia> extrema(mapslices(std, z; dims = 2))
(0.999483368783979, 1.0003834076282714)

julia> quantile(z_dens, ps)
11-element Array{Float64,1}:
 -0.00012070140687908017
 -0.0001063122814842569 
 -0.00010400024477877038
 -0.00010240952128016883
 -0.00010104242015569829
 -9.979156638516813e-5  
 -9.854879329589064e-5  
 -9.727505743865635e-5  
 -9.584836020836684e-5  
 -9.393933488262388e-5  
 -8.430228880599951e-5  

julia> quantile(sim_z_dens, ps)
11-element Array{Float64,1}:
 -0.00012556606255873136
 -0.00010655950256849811
 -0.0001041351251443149 
 -0.0001024130143312689 
 -0.00010104015366172799
 -9.97620830194371e-5   
 -9.849225327092044e-5  
 -9.72325802687073e-5   
 -9.577554898483448e-5  
 -9.377537093476237e-5  
 -8.342424675283704e-5  

```

---

<div class="post-metadata">

**Author:** ![jocklawrie](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jocklawrie/32/33616_2.png) [@jocklawrie](https://discourse.julialang.org/u/jocklawrie)\
**Post date:** [March 7, 2019, 4:54am UTC](https://discourse.julialang.org/t/generating-high-dimensional-sobol-sequences/21521/3 "2019-03-07T04:54:42Z")

</div>

Thanks for responding.

I don’t think numerical error is at play here. While it is true that across each of the 100 dimensions we get a good (low-discrepancy) sample from the univariate standard normal, I think we are seeing a natural consequence of high-dimensional sampling…that is, a handful of points dominating the distribution in terms of density.

I suspect this is unavoidable because within a draw all of the 100 component densities are less than 0.4 (`pdf(Normal(0, 1), 0)`), and their product is usually very small. But sometimes (rarely) most of the 100 component densities will be relatively close to 0.4 (or at least on the order of 10^-1), which results in a much larger density than that of most draws…several orders of magnitude in fact.