Vtyjkyuk

Suppose we have several features (e.g. $geq20$) that do not follow a Gaussian distribution. Do we have to worry about the features not following a Gaussian distribution if we apply standardization on the data?

Namely, even if the features do not follow a normal distribution initially, aren't they made to follow Gaussian distribution after standardization with mean $0$ and variance $1$?

The short answer: yes, you do need to worry about your data's distribution not being normal, because standardization does not transform the underlying distribution structure of the data. If $XsimmathcalN(mu, sigma^2)$ then you can transform this to a standard normal by standardizing: $Y:=(X-mu)/sigma simmathcalN(0,1)$. However, this is possible because $X$ already follows a normal distribution in the first place. If $X$ has a distribution other than normal, standardizing it in the same way as above will generally not make the data normally distributed.

A simple example of exponentially distributed data and its standardized version:

x <- rexp(5000, rate = 0.5)
y <- (x-mean(x))/sd(x)
par(mfrow = c(2,1))
hist(x, freq = FALSE, col = "blue", breaks = 100, xlim = c(min(x), quantile(x, 0.995)),
 main = "Histogram of exponentially distributed data X with rate = 0.5")
hist(y, freq = FALSE, col = "yellow", breaks = 100, xlim = c(min(y), quantile(y, 0.995)),
 main = "Histogram of standardized data Y = ( X-E(X) ) / StDev(X)")

Now if we check the mean and standard deviation of the original data $x$, we get

c(mean(x), sd(x))
[1] 2.044074 2.051816

whereas for the standardized data $y$, the corresponding results are

c(mean(y), sd(y))
[1] 7.136221e-17 1.000000

As you can see, the distribution of the data after standardization is decidedly not normal, even though the mean is (practically) 0 and the variance 1. In other words, if the features do not follow a normal distribution before standardization, they will not follow it after the standardization either.