Why does Central Limit Theorem break down in my simulation? Let say I have following numbers:
4,3,5,6,5,3,4,2,5,4,3,6,5

I sample some of them, say, 5 of them, and calculate the sum of 5 samples.
Then I repeat that over and over to get many sums, and I plot the values of sums in a histogram, which will be Gaussian due to the Central Limit Theorem.
But when they are following numbers, I just replaced 4 with some big number:
4,3,5,6,5,3,10000000,2,5,4,3,6,5

Sampling sums of 5 samples from these never becomes Gaussian in histogram, but more like a split and becomes two Gaussians. Why is that?
 A: I'd just like to explain, using complex cumulant-generating functions, why everyone keeps blaming this on skew.
Let's write the random variable you're sampling as $\mu+\sigma Z$, where $\mu$ is the mean and $\sigma$ the standard deviation so $Z$ has mean $0$ and variance $1$. The cumulant-generating function of $Z$ is $-\frac{1}{2}t^2-\frac{i\gamma_1}{6}t^3+o(t^3)$. Here $\gamma_1$ denotes the skew of $Z$; we could write it in terms of the skew $\kappa_3$ of the original variable $\mu+\sigma Z$, viz. $\gamma_1=\sigma^{-3}\kappa_3$. 
If we divide the sum of $n$ samples of $Z$'s distribution by $\sqrt{n}$, the result has cgf $$n\left(-\frac{1}{2}\left(\frac{t}{\sqrt{n}}\right)^2-\frac{i\gamma_1}{6}\left(\frac{t}{\sqrt{n}}\right)^3\right)+o(t^3)=-\frac{1}{2}t^2-\frac{i\gamma_1}{6\sqrt{n}}t^3+o(t^3).$$For a Normal approximation to be valid at large enough $t$ for the graph to look right, we need sufficiently large $n$. This calculation motivates $n\propto\gamma_1^2$. The two samples you considered have very different values of $\gamma_1$.
A: Let's recall, precisely, what the central limit theorem says.

If $X_1, X_2, \cdots, X_k$ are independent and identically distributed random variables with (shared) mean $\mu$ and standard deviation $\sigma$, then $\frac{X_1 + X_2 + \cdots + X_k}{k\frac{\sigma}{\sqrt{k}}}$ converges in distribution to a standard normal distribution $N(0, 1)$ (*).

This is often used in the "informal" form:

If $X_1, X_2, \cdots, X_k$ are independent and identically distributed random variables with (shared) mean $\mu$ and standard deviation $\sigma$, then $X_1 + X_2 + \cdots + X_k$ converges "in distribution" to a standard normal distribution $N(k \mu, \sqrt{k} \sigma)$.

There's no good way to make that form of the CLT mathematically precise, since the "limit" distribution change, but it's useful in practices.
When we have a static list of numbers like
4,3,5,6,5,3,10000000,2,5,4,3,6,5

and we are sampling by taking a number at random from this list, to apply the central limit theorem we need to be sure that our sampling scheme satisfies these two conditions of independence and identically distributed.


*

*Identically distributed is no problem: each number in the list is equally likely to be chosen.

*Independent is more subtle, and depends on our sampling scheme.  If we are sampling without replacement, then we violate independence.  It is only when we sample with replacement that the central limit theorem is applicable.


So, if we use with replacement sampling in your scheme, then we should be able to apply the central limit theorem.  At the same time, you are right, if our sample is of size 5, then we are going to see very different behaviour depending on if the very large number is chosen, or not chosen in our sample.
So what's the rub?  Well, the rate of convergence to a normal distribution is very dependent on the shape of the population we are sampling from, in particular, if our population is very skew, we expect it to take a long time to converge to the normal.  This is the case in our example, so we should not expect that a sample of size 5 is sufficient to show the normal structure.

Above I repeated your experiment (with replacement sampling) for samples of size 5, 100, and 1000.  You can see that the normal structure is emergent for very large samples. 
(*) Note there are some technical conditions needed here, like finite mean and variance.  They are easily verified to be true in our sampling from a list example.
A: In general, the size of each sample should be more than $5$ for the CLT approximation to be good. A rule of thumb is a sample of size $30$ or more. But, with the population of your first example, $5$ is OK.
pop <- c(4, 3, 5, 6, 5, 3, 4, 2, 5, 4, 3, 6, 5)
N <- 10^5
n <- 5
x <- matrix(sample(pop, size = N*n, replace = TRUE), nrow = N)
x_bar <- rowMeans(x)
hist(x_bar, freq = FALSE, col = "cyan")
f <- function(t) dnorm(t, mean = mean(pop), sd = sd(pop)/sqrt(n))
curve(f, add = TRUE, lwd = 2, col = "red")


In your second example, because of the shape of the population distribution (for one thing, it's too much skewed; read the comments by guy and Glen_b bellow), even samples of size $30$ won't give you a good approximation for the distribution of the sample mean using the CLT.
pop <- c(4, 3, 5, 6, 5, 3, 10000000, 2, 5, 4, 3, 6, 5)
N <- 10^5
n <- 30
x <- matrix(sample(pop, size = N*n, replace = TRUE), nrow = N)
x_bar <- rowMeans(x)
hist(x_bar, freq = FALSE, col = "cyan")
f <- function(t) dnorm(t, mean = mean(pop), sd = sd(pop)/sqrt(n))
curve(f, add = TRUE, lwd = 2, col = "red")


But, with this second population, samples of, say, size $100$ are fine.
pop <- c(4, 3, 5, 6, 5, 3, 10000000, 2, 5, 4, 3, 6, 5)
N <- 10^5
n <- 100
x <- matrix(sample(pop, size = N*n, replace = TRUE), nrow = N)
x_bar <- rowMeans(x)
hist(x_bar, freq = FALSE, col = "cyan")
f <- function(t) dnorm(t, mean = mean(pop), sd = sd(pop)/sqrt(n))
curve(f, add = TRUE, lwd = 2, col = "red")


A: Short answer is, you don't have a big enough sample to make central limit theorem apply.
