Calculating the likelihood of a whole dataset Of my dataset $X$ I want to calculate the likelihood using my GMM model (for BIC). Mathematically it seems to make sense that as the samples are independant $P(X) = P(X_1)P(X_2)..$, so I would get the likelihood by taking the product of the likelihoods per sample.
But it doesn't seem intuitively correct that the more samples I have the less likely my model. It also means that if just one sample that doesn't fit will crash my likelihood, it therefore seems to make more sense (to me) to average the likelihoods per sample.
What's the right thing to do here (I feel I'm messing some very basic statistics up)?
 A: 
But it doesn't seem intuitively correct that the more samples I have the less likely my model.

I think your intuition has it the wrong way around. You're not computing how likely your model is - you're computing how likely the data are, GIVEN your model. You do this by multiplying together $n$ probabilities in the unit interval, so you should expect likelihood to decrease exponentially with $n$.
The key insight is that you don't use measures like BIC or AIC as a one to one measure of how likely your model is to be true. You compare different models by comparing their scores. For instance, the difference between BIC scores of models A and B, is an approximation of the ratio of the probabilities of the models being true (assuming the true model is in the set of candidate models). Hence, it doesn't really matter that the likelihood of the true model decreases with n, as long as the likelihoods for other models drop quicker, so the difference between their BICs grow.
A: 
it doesn't seem intuitively correct that the more samples I have the less likely my model

Maximum likelihood is used to find parameters of the model, so it isn't really a problem.

seems to make more sense (to me) to average the likelihoods per sample.

You're right if you would want to make likelihood more meaningful.
But it doesn't really matter for Maximum Likelihood since it is used for finding $$\hat{\theta} = \underset{\theta}{argmax} \underset{i < N}{\prod}{P(x_i; \theta)}$$  $$=\underset{\theta}{argmax} \sqrt[n]{\underset{i < N}{\prod}{P(x_i; \theta)}}$$ (nth root is increasing function, and correct averaging for product would be geometric averaging.).

It also means that if just one sample that doesn't fit will crash my likelihood

What do you mean by 'crashing likelihood'? If you mean that your likelihood becomes zero, this is impossible for GMM, since gaussian mixture's PDF is nonzero everywhere.
