function x = chi2_inv(p,k)
% x = chi2_inv(F,k)
% compute a chi-squared statistic from a distribution with k degrees of freedom
% and a probability level of p
% using the Wilson-Hilferty transformation
% https://en.wikipedia.org/wiki/Chi-squared_distribution#Asymptotic_properties 

% Wilson-Hilferty transformation

m = 1 - 2 / (9*k);    % mean of cube root of X / k
s = sqrt(2/(9*k));    % standard deviation of cube root of X / k

z = norm_inv ( p, m, s);

x = k * z.^3;

% ----------------------------------------------------------------------
