FSolve
Description
FSolve solves for the coefficients of Fleishman’s power method, which turns normal random numbers into ones with a requested skewness and excess kurtosis. It uses Newton’s method (up to 50 iterations) and returns the three coefficients, the number of iterations and the residual, so the caller can tell whether it converged. It doesn’t check whether the requested skewness and kurtosis are possible; see Fkmin.
Used by: MCM(_V), MCM_CI(_V), MCM_error(_V), rerror, rerrorV
Source Code
FSolve = LAMBDA(
skew, // requested skewness
kurt, // requested excess kurtosis (not checked against min required)
[b], [c], [d], [n], // for recursion
LET(
b0, IF(ISOMITTED(b), 1, b),
c0, IF(ISOMITTED(c), skew/6, c),
d0, IF(ISOMITTED(d), 0, d),
n0, IF(ISOMITTED(n), 0, n),
F_1, b0^2+6*b0*d0+2*c0^2+15*d0^2-1,
F_2, 2*c0*(b0^2+24*b0*d0+105*d0^2+2)-skew,
F_3, 24*(b0*d0+c0^2*(1+b0^2+28*b0*d0)+d0^2*(12+48*b0*d0+141*c0^2+225*d0^2))-kurt,
resid, SQRT(F_1^2+F_2^2+F_3^2),
IF(OR(n0 >= 50, resid < 0.0000001),
HSTACK(b0, c0, d0, n0, resid),
LET(
J, VSTACK(
HSTACK(2*b0+6*d0, 4*c0, 6*b0+30*d0),
HSTACK(4*c0*b0+48*c0*d0, 2*b0^2+48*b0*d0+210*d0^2+4, 48*c0*b0+420*c0*d0),
HSTACK(24*(d0+2*b0*c0^2+28*d0*c0^2+48*d0^3),
24*(2*c0+2*b0^2*c0+56*b0*d0*c0+282*c0*d0^2),
24*(b0+28*b0*c0^2+24*d0+144*b0*d0^2+282*c0^2*d0+900*d0^3))
),
Fvec, VSTACK(F_1,F_2,F_3),
delta, MMULT(MINVERSE(J), Fvec),
FSolve(skew, kurt, b0-INDEX(delta,1,1), c0-INDEX(delta,2,1), d0-INDEX(delta,3,1), n0+1)
)
)
)
);