Cannot figure out simple use of Cumulants.jl

Question

I cannot for the life of me figure out how to use Cumulants.jl to get moments or cumulants from some data. I find the docs (https://juliahub.com/docs/Cumulants/Vrq25/1.0.4/) completely over my head.

Suppose I have a vector of some data e.g.:

using Distributions
d = rand(Exponential(1), 1000)

The documentation suggests, so far as I can understand it, that cumulants(d, 3) should return the first three cumulants. The function is defined like so:

cumulants(data::Matrix{T}, m::Int = 4, b::Int = 2) where T<: AbstractFloat

a Matrix in Julia is, so far as I understand, a 2D array. So I convert my data to a 2D array:

dm = reshape(d, length(d), 1)

But I get:

julia> cumulants(dm,3)
ERROR: DimensionMismatch("bad block size 2 > 1")

My question concisely: how do I use Cumulants.jl to get the first m cumulants and the first m moments from some simulated data?

Thanks!

EDIT: In the above example, c = cumulants(dm,3,1) as suggested in a comment will give, for c:

3-element Array{SymmetricTensors.SymmetricTensor{Float64,N} where N,1}:
 SymmetricTensors.SymmetricTensor{Float64,1}(Union{Nothing, Array{Float64,1}}[[1.0122452678071678]], 1, 1, 1, true)
 SymmetricTensors.SymmetricTensor{Float64,2}(Union{Nothing, Array{Float64,2}}[[1.0336298356976195]], 1, 1, 1, true)
 SymmetricTensors.SymmetricTensor{Float64,3}(Union{Nothing, Array{Float64,3}}[[2.5438037582591146]], 1, 1, 1, true)

I find that I can access the first, second, and third cumulants by:

c[1][1]
c[2][1,1]
c[3][1,1,1]

Which I arrived at essentially by guessing. I have no idea why this nutty output format exists. I still cannot figure out how to get the first m cumulants as a vector easily.

Answer 1

As I wrote in the comments, if you have a univariate problem you should use cumulants(dm,3,1) as the cumulants are calulated using tensors and the tensors are saved in a block structure, where the blocks are of size bxb, i.e. the third argument in the function call. However, If you have only one column, the size of the tensors will be 1, so that it doesn't make sense to save it in a 2x2 block.

To access the cumulants in Array form you have to convert them first. This is done by Array(cumulant(data, nc, b)[c]), where nc is the number of cumulants you want to calculate, b is the block size (for efficient storage of the tensors), and c is the cumulant you need. Summing up:

using Cumulants

# univariate data

unidata = rand(1000,1)
uc = cumulants(unidata, 3, 1)
Array(uc[1])
#1-element Array{Float64,1}:
# 0.48772026299259374
Array(uc[2])
#1×1 Array{Float64,2}:
# 0.0811428357438324
Array(uc[3])
#[:, :, 1] =
# 0.0008653019738796724

# multivariate data

multidata = rand(1000,3)
mc = cumulants(multidata, 3, 2)
Array(mc[1])
#3-element Array{Float64,1}:
# 0.5024511157116442
# 0.4904838734508787
# 0.48286680648519215
Array(mc[2])
#3×3 Array{Float64,2}:
#  0.0834021   -0.00368562  -0.00151614
# -0.00368562   0.0835084    0.00233202
# -0.00151614   0.00233202   0.0808521
Array(mc[3])
# [:, :, 1] =
#  -0.000506926  -0.000763061  -0.00183751
#  -0.000763061  -0.00104804   -0.00117227
#  -0.00183751   -0.00117227    0.00112968
# 
# [:, :, 2] =
#  -0.000763061  -0.00104804   -0.00117227
#  -0.00104804    0.000889305  -0.00116559
#  -0.00117227   -0.00116559   -0.000106866
# 
# [:, :, 3] =
#  -0.00183751  -0.00117227    0.00112968
#  -0.00117227  -0.00116559   -0.000106866
#   0.00112968  -0.000106866   0.00131965

The optimal size of the blocks can be found in their software paper (https://arxiv.org/pdf/1701.05420.pdf), where they write (for proper latex formatting have a look at the paper):

5.2.1. The optimal size of blocks. The number of coefficients required to store a super-symmetric tensor of order d and n dimensions is equal to (d+n−1 over n). The storage of tensor disregarding the super-symmetry requires n^d coefficients. The block structure introduced in [49] uses more than minimal amount of memory but allows for easier further processing of super-symmetric tensors.If we store the super-symmetric tensor in the block structure, the block size parameter b appears. In our implementation in order to store a super-symmetric tensor in the block structure we need, assuming n|b, an array of (n over b)^d pointers to blocks and an array of the same size of flags that contain the information if a pointer points to a valid block. Recall that diagonal blocks contain redundant information.Therefore on the one hand, the smaller the value of b, the less redundant elements on diagonals of the block structure. On the other hand, the larger the value of b,the smaller the number of blocks, the smaller the blocks’ operation overhead, and the fewer the number of pointers pointing to empty blocks. For detailed discussion of memory usage see [49]. The analysis of the influence of the parameter b on the computational time of cumulants for some parameters are presented in Fig. 2. We obtain the shortest computation time for b = 2 in almost all test cases, and this value will be set as default and used in all efficiency tests. Note that for b = 1we loose all the memory savings.

Answer 2

Using Oskar's helpful answer, I thought I'd provide my wrapper function which accomplishes the goal of returning a vector of the first m cumulants, given an input of a 1D array of data.

using Cumulants
function mycumulants(d, m) # given a 1D array of data d, return a vector of the first m cumulants
    res = zeros(m)
    dm = reshape(d, length(d), 1) # Convert 1D array to 2D
    c = cumulants(dm, m, 1) # Need the 1 (block size) or else it errors
    for i in 1:m
        res[i] = Array(c[i])[1]
    end
    return(res)
end

But it turns out this is really really slow compared to just directly calculating raw moments and coverting them to cumulants by e.g. k[5] = u[5] - 5*u[4]*u[1] - 10*u[3]*u[2] + 20*u[3]*u[1]^2 + 30*u[2]^2*u[1] - 60*u[2]*u[1]^3 + 24*u[1]^5 so I think I won't be using Cumulants.jl after all for my purposes, which only involve univariate data at this time.

Example of time difference for calculating the first six cumulants from some simulated data:

----Data set 2----

Direct calculation:
  1.997 ms (14 allocations: 469.47 KiB)

Cumulants.jl:
  152.798 ms (318435 allocations: 17.59 MiB)