# [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 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] 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[i] * prod([S[j] for j in repr[i]]) for i in 1:size(repr,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 == S/N
@assert result == (N * S - S^2)/(N*(N-1))
@assert result == (2S^3 - 3N*S*S +N^2 * S) / (N*(N-1)*(N-2))
@assert result == (-6*S^4 + 12 * N * S^2 * S - 3 * N * (N-1) * S^2 - 4 * N * (N+1) * S * S + N^2*(N+1)*S)/(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