Sorry, I should have been clear that the toolbox was only compatible with ArcGIS 10. You are on the right track with calculate statistics. Make sure you set the skip value to 1, then add the raster to ArcMap (or navigate to it in ArcCatalog) and right click on it and select properties. In the layer properties box select the source tab and scroll down to the bottom and you will see the statistics section and can retrieve the min, max, mean and stdv (can copy and paste the values). You can then use the values to perform a normalization on your raster in the raster calculator (Arctoolbox > Spatial Analyst Tools > Map Algebra > Raster Calculator) with the following syntax.
(x - mean) / stdv
where; x is your raster, mean and stdv are the real values of the respective statistical moments.
The above formula from your original post does not transform to a standard variable space. This transformation is intended to scale the mean to 0 and stdv to 1 while maintaining the shape of the original distribution. The resulting data range is dictated by the range in your original data. This formula, in effect, makes the negative and positive bounds symmetrical (e.g., -1 to 1). If the original data distribution is non-normal the results can be unexpected. If you just want your data in the same scale and it is all positive, you could just perform a "row standardization" by dividing your raster by its max value to (mostly) return a range of 0-1. Here is a method that accounts for negative values and also reliably returns a range of 0-1. ( x - min(x) ) / ( max(x) - min(x) )
If you have access to R, the behavior of these three transformations can be readily observed given a normal and skewed distribution. Note that you can plot the distribution of any of the resulting transformations and it does not change shape, just range. Here is the code (just copy and paste).
############################
# Based on a non-normal distribution
############################
x <-runif(100,1,100)
summary(x)
plot(density(x))
x1 <- ( x - mean(x) ) / sd(x)
summary(x1)
x2 <- x / max(x)
summary(x2)
x3 <- ( x - min(x) ) / ( max(x) - min(x) )
summary(x3)
# When distribution has negative values (note that the regular row standardization [x/max(x)] does not scale correctly)
x[1] <- -10
x2 <- x / max(x)
x3 <- ( x - min(x) ) / ( max(x) - min(x) )
summary(x2);summary(x3)
############################
# Based on a non-normal distribution
############################
x <- rweibull(1e5,1.5,33)
summary(x)
plot(density(x))
x1 <- ( x - mean(x)) / sd(x)
summary(x1)
x2 <- x / max(x)
summary(x2)
x3 <- ( x - min(x) ) / ( max(x) - min(x) )
summary(x3)