[Help needed] Translate 10 lines from Maple to Julia

Hey,

I have a somewhat weird request: I do not speak Maple at all, but i need to translate the 10 lines of code from this document, to be able to use the function ks() that they created on page 2. Notes on underlying other functions are on page 3 and 4, the rest of the document gives generalisations that I do not need.

This function computes a polynomial in S_1,...,S_m that corresponds to the expression of the k-statistic (unbiaised estimator of the kth cumulant):

Details on the corresponding theory, at least for the first few cases, is available there: k-Statistic -- from Wolfram MathWorld

By reading the Maple code, do you think it’s going to be complicated ?

Edit: I found the partition and stirling2 functions in the Combinatorics parckage. This is what i have for the moment:

import Combinatorics


# The partition Maple function is replaed by : 
# Combinatorics.integer_partitions(k)

# stirling2 is replaced by 
# Combinatorics.stirlings2(k,n)


# Let's try maketab: 
function makeTab_ss(N)
    

end

I still need to undertand wtf is doing the makeTab_ss function :wink:

1 Like

If you need just naive Julia implementation (not-as-fast-as-it-could-be) only problem is to decipher Maple code. I quite rusty in things like Maple, but, if my time allows, I can try to help you with that. You can send me message if you want so more help.

I need something that works and i do not need it to be fast for the moment, as it would be a one time computation that I then store somewhere.

If you can help me understand the Maple code it would be great !

Ok, how much time we have? I’m now quite busy, but I probably can sit down in few hours. Here is now 16:08, around 19 I should have some time.

Well, i’ll be happy to show you what i have by 19 ! Thanks a lot.

Ok i do have the following code:

# The goal of this file is to compute the k_statistics polynomials ad given by 
# https://discourse.julialang.org/t/help-needed-translate-10-lines-from-maple-to-julia/60650
# and 
# https://iris.unito.it/retrieve/handle/2318/1561388/399488/FASTpolykaysMAPLE.pdf
# also see 
# Di Nardo E., G. Guarino, D. Senato (2007), A new method for fast computing unbiased estimators of cumulants. In press Statistics and Computing. http://www.springer.com/statistics/computational/journal/11222 (download from http://www.unibas.it/utenti/dinardo/lavori.html)


import Combinatorics
import Polynomials
import StatsBase
import Random
import TypedPolynomials

function makeTab_ss(N::T) where T<:Integer
    part = Combinatorics.integer_partitions(N) 
    dict = Dict{Int64, Int64}()
    Nf = factorial(N)
    result = Array{Tuple{Rational{Int64}, Vector{Int64}}}(undef,0)
    for y in part
        denom = T(1)
        dict = StatsBase.countmap(y)
        for (val,numoccur) in dict
            denom *= factorial(val)^numoccur * factorial(numoccur)
        end
        push!(result,(Nf//denom,y))
    end
    return result
end

function Stirling_triangle(N::T) where T<:Integer
    arr = zeros(T,N,N)
    for i in 1:N
        for j in 1:i
            arr[i,j] = Combinatorics.stirlings2(i,j)
        end
    end
    return arr
end

function makeK_s(N::T) where T<:Integer
    result = Array{Polynomials.Polynomial{Float64}}(undef,0)
    Stirling2 = Stirling_triangle(N)
    for i in 1:N
        rez = zeros(i+1)
        for j in 1:i
            rez[j+1] = Combinatorics.stirlings2(i,j) * (-1)^(j-1) * factorial(j-1)
        end
        push!(result,Polynomials.Polynomial(rez))
    end
    return result
end


function expand_eval(v,u)
    result = Array{eltype(u)}(undef,0)
    for (number,array) in v
        rez = number * prod(u[i] for i in array)::eltype(u)
        push!(result,rez)
    end
    return result
end

function fd(j,k)
    Polynomials.fromroots([i+k for i in 0:(j-k)])
end
# checks from the note in the source: 
@assert fd(3,0) == Polynomials.Polynomial([0,-6,11,-6,1])
@assert fd(3,1) == Polynomials.Polynomial([-6,11,-6,1])
@assert fd(3,2) == Polynomials.Polynomial([6,-5,1])

function build_u2(N)
    result = []
    for i in 1:N
        rez = (-1)^(i-1)*factorial(i-1)*fd(N-1,i)
        push!(result,rez)
    end
    return result
end

function compute_k_representation(n::T,k_idx) where T

    v = makeTab_ss(k_idx)
    u = makeK_s(k_idx)
    v2 = expand_eval(v,u)
    poly = zero(eltype(u))
    u2 = build_u2(k_idx)
    common_denominator = prod(n - x for x in 0:(k_idx-1))
    factors = []
    for j in 1:size(v2,1)
        poly = sum((v2[j].coeffs[i+1]*u2[i]) for i in 1:(length(v2[j].coeffs)-1))
        numerator = poly(n)
        push!(factors,numerator/common_denominator)
    end
    which_S = [v[i][2] for i in 1:length(v)]
    return (factors,which_S)
end

function compute_k(N,S)
    m = length(S)
    result = []
    for k in 1:m
        repr = compute_k_representation(N,k)
        rez = sum(repr[1][i] * prod([S[j] for j in repr[2][i]]) for i in 1:size(repr[1],1))
        push!(result,rez)
    end
    return result
end

# Simple experiment: 
# Vary the number of observations N, and the mean and variance of the normal random vector. 
# The first cumulant should be the mean, 
# the second should be the variance 
# and all others should be 0 since it's a normal. 
# the approximation should get better and better as you increase the number of pointS. 
mean = 0
var = 1
m=7
N_obs = 100000
Random.seed!(123)
data = mean .+ sqrt(var) .* randn(N_obs)
S_obs = [sum(data .^ i) for i in 1:m]
cumulants_direct = compute_k(big.(N_obs),big.(S_obs)) # directly without the poynomials.

# As a simple check, we can verify that the polynomials holds : 
TypedPolynomials.@polyvar N
TypedPolynomials.@polyvar S[1:7] ##########" put m here in a fixed way, as the infrastructure does not allow for variables.
result = compute_k(N,S);
@assert result[1] == S[1]/N
@assert result[2] == (N * S[2] - S[1]^2)/(N*(N-1))
@assert result[3] == (2S[1]^3 - 3N*S[1]*S[2] +N^2 * S[3]) / (N*(N-1)*(N-2))
@assert result[4] == (-6*S[1]^4 + 12 * N * S[1]^2 * S[2] - 3 * N * (N-1) * S[2]^2 - 4 * N * (N+1) * S[1] * S[3] + N^2*(N+1)*S[4])/(N*(N-1)*(N-2)*(N-3))
cumulants_poly = [result[i](N=>big.(N_obs),S=>big.(S_obs)) for i in 1:m]


display([cumulants_direct cumulants_poly])

The code is clearly suboptimal, however with the introduction of TypedPolynomials, i was able to check that it is correct.

Now the question is :

  • Is it possible to compute once and for all the polynomials, for any N and S, and store them to disk somewhere in a package that i can load easily ?
  • Is it possible to do the same thing (pre-compute the polynomials) for a given N, on the fly ?

The use case is that i might compute the value of these polynomials a lot of time, for varying S and for fixed N.

1 Like