Computing the real solutions of Fleishman's equations for simulating non-normal data

Fleishman's power method is frequently used to simulate non-normal data with a desired skewness and kurtosis. Fleishman's method requires solving a system of nonlinear equations to find the third-order polynomial weights that transform a standard normal variable into a non-normal variable with desired moments. Most users of the power method seem unaware that Fleishman's equations have multiple solutions for typical combinations of skewness and kurtosis. Furthermore, researchers lack a simple method for exploring the multiple solutions of Fleishman's equations, so most applications only consider a single solution. In this paper, we propose novel methods for finding all real-valued solutions of Fleishman's equations. Additionally, we characterize the solutions in terms of differences in higher order moments. Our theoretical analysis of the power method reveals that there typically exists two solutions of Fleishman's equations that have noteworthy differences in higher order moments. Using simulated examples, we demonstrate that these differences can have remarkable effects on the shape of the non-normal distribution, as well as the sampling distributions of statistics calculated from the data. Some considerations for choosing a solution are discussed, and some recommendations for improved reporting standards are provided.