We introduce a new approximation method for the distribution of functions of random variables that are real-valued. The approximation involves moment matching and exploits properties of the class of normal inverse Gaussian distributions. In the paper we examine the how well the different approximation methods can capture the tail behavior of a function of random variables relative each other. This is done by simulate a number functions of random variables and then investigate the tail behavior for each method. Further we also focus on the regions of unimodality and positive definiteness of the different approximation methods. We show that the new method provides equal or better approximations than Gram-Charlier and Edgeworth expansions.