How to transform data to normality? We have financial some data (500-1000 samples), which is not normally distributed (well known fact from the literature). I have some ideas to do parametric transformations of this data (using some other data) to produce "adjusted" series. My goal is to find a transformation that makes the series normally distributed (with mean 0 and std deviation 1). What is the most appropriate statistic and corresponding test to optimize my parameters and determine if the outcome can be considered normally distributed?
Please also point me to an implementation, ideally in C/C++ or java.
 A: For financial data I have successfully used heavy-tail Lambert W x Gaussian transformations.


*

*Pyhon: gaussianize is an sklearn-type implementation of the IGMM algorithm in Python.

*C++: the lamW R package has an elegant (and fast) C++ implementation of Lambert's W function.  This can be a starting point for a full C++ implementation of IGMM or MLE for Lambert W x Gaussian transformations.

*R: the LambertW R package is a full implementation of the Lambert W x F framework (simulation, estimation, plotting, transformation, testing).
As an illustration consider the SP500 return series in R.
library(MASS)
data(SP500)
yy <- ts(SP500)
library(LambertW)
test_norm(yy)


## $seed
## [1] 516797
## 
## $shapiro.wilk
## 
##  Shapiro-Wilk normality test
## 
## data:  data.test
## W = 1, p-value <2e-16
## 
## 
## $shapiro.francia
## 
##  Shapiro-Francia normality test
## 
## data:  data.test
## W = 1, p-value <2e-16
## 
## 
## $anderson.darling
## 
##  Anderson-Darling normality test
## 
## data:  data
## A = 20, p-value <2e-16

As is well-known, financial data typically have fat tails and are sometimes negatively skewed.  For the SP500 case, skewness is not too large, but it exhibits high kurtosis (7.7).  Also several normality tests clearly reject the null hypothesis of a marginal Gaussian distribution.
Since we only have to deal with heavy tails, but not skewness, let's fit a heavy-tailed Lambert W x Gaussian distribution using a method of moments estimator (one could also use the maximum likelihood estimator (MLE) with MLE_LambertW()).
# fit a heavy tailed Lambert W x Gaussian
mod <- IGMM(yy, type = "h")
mod

## Call: IGMM(y = yy, type = "h")
## Estimation method:  IGMM 
## Input distribution:  Any distribution with finite mean & variance and kurtosis = 3. 
## mean-variance Lambert W x F type ('h' same tails; 'hh' different tails; 's' skewed):  h 
## 
##  Parameter estimates:
##    mu_x sigma_x   delta 
##    0.05    0.72    0.16 
## 
##  Obtained after 4 iterations.

The heavy tail parameter $\widehat{\delta} = 0.16$ is significantly different from zero and implies heavy tails.  For such a $\delta$ moments up to order $1 / \widehat{\delta} = 6.29$ exist. 
The model check question is of course if the back-transformed data does indeed have a Gaussian distribution.  Let's check again using test_norm():
# transform data to input data (which presumably should have Normal distribution); use return.u = TRUE to get zero-mean, unit variance data
xx <- get_input(mod, return.u = FALSE)
test_norm(xx)


## $seed
## [1] 268951
## 
## $shapiro.wilk
## 
##  Shapiro-Wilk normality test
## 
## data:  data.test
## W = 1, p-value = 0.2
## 
## 
## $shapiro.francia
## 
##  Shapiro-Francia normality test
## 
## data:  data.test
## W = 1, p-value = 0.2
## 
## 
## $anderson.darling
## 
##  Anderson-Darling normality test
## 
## data:  data
## A = 0.7, p-value = 0.07

I think the plot and normality test results speak for themselves.
The package also provides a single function that does all these steps at once: Gaussianize() (this is also what the Python package implements).
A: It appears that you are just asking for a test for normality.  If so, Shapiro-Wilk is hard to beat.  This is not, however, the easiest test in the pantheon to implement.  
Why not just use R?  The shapiro.test function will do the work for you.  
