How would I go about calculating the probability that the observed mean of a collection of data generated by a random process, which is governed by a particular statistical model, would deviate from the theoretical mean by a given amount?
For example, the sum of the rolls of two six-sided dice in theory follows a statistical model that resembles a Gaussian probability distribution curve, between the numbers 2 and 12, with a mean of 7. If I was to roll and sum the results of two dice 100 times, how can I calculate the probability that the observed mean of data would be 6, rather than 7?
In other words, if I was to do that test 100 times and the observed mean was calculated to be 6, how can I calculate the probability that the deviation in the mean could be because of random chance, as opposed to the underlying (assumed) statistical model being incorrect (i.e. the dice being biased)?
I know some basic statistics, but I'm far from an expert. I would think that the observed mean of a particular statistical sample would have a probability distribution of its own, but I'm not sure how it could be calculated.
Let the random variable $X_i$ denote the outcome of the $i$-th dice roll. Then the observed mean $\bar X$ of $n$ dice rolls will simply be their sum $\Sigma_X$ divided by the number of rolls $n$, i.e. $$\bar X = \frac{\Sigma_X}n = \frac{X_1 + X_2 + \dots + X_n}n = \frac1n \sum_{i=1}^n X_i.$$
This implies that the distribution of $\bar X$ is simply the distribution of $\Sigma_X$ transformed by uniformly scaling the values of the result by $\frac1n$. In particular, the probability that $\bar X \le 6$, for example, is equal to the probability that $\Sigma_X \le 6n$.
Now, in your case, each $X_i$ is the sum of the rolls of two six-sided dice, so $\Sigma_X$ is simply the sum of $2n$ six-sided dice. This distribution does not (as far as I know) have a common name or any particularly nice algebraic form, although it is closely related to a multinomial distribution with $2n$ trials and six equally likely outcomes per trial. However, it's not particularly hard to calculate it numerically, e.g. by starting with the singular distribution with all probability mass concentrated at zero, and convolving it $2n$ times with the uniform distribution over $\{1,2,3,4,5,6\}$.
For example, this Python program calculates and prints out the distribution of the sum of 200 six-sided dice, both as the raw point probability of rolling each sum and as the cumulative probability of rolling each sum or less. From the line numbered 600, we can read that the probability of rolling a sum of 600 or less (and thus a mean of 6 or less over 100 pairs of rolls) is approximately 0.00001771.
(Alternatively, we could use a tool made specifically for this job, although we'd need to switch to the export view to get an accurate numerical value for a result so far away from the expected value.)
For large numbers of rolls (say, more than 10) we can also get pretty good results simply by approximating the distribution of $\bar X$ (and/or $\Sigma_X$) with a normal distribution with the same expectation and variance. This is a consequence of the central limit theorem, which says that the distribution of the sum (and the average) of a large number of identically distributed random variables tends towards a normal distribution.
By the linearity of the expected value, the expected value of $\Sigma_X$ is simply the sum of the expected values of the dice rolls, and the expected value of $\bar X$ is their average (i.e. the same as for each individual roll, if the rolls are identically distributed). Further, assuming that the dice rolls are independent (or at least uncorrelated), the variance of $\Sigma_X$ will be the sum of the variances of the dice rolls, and the variance of $\bar X$ will be $\frac1{n^2}$ of the variance of $\Sigma_X$ (i.e. respectively $n$ and $\frac1n$ times the variance of each individual roll, if they're all the same).
Since the variance of a single $k$-sided die roll (i.e. a uniformly distributed random variable over $\{1, 2, \dots, k\}$ is $(k^2-1)\mathbin/12$, and the variance of the sum of two rolls is simply twice that, it follows that the variance of a roll of two six-sided dice is $2 \times 35 \mathbin/ 12 = 35 \mathbin/ 6 \approx 5.83$. Therefore, over 100 such double rolls, the distribution of the sum $\Sigma_X$ is well approximated by a normal distribution with mean $\mu = 700$ and variance $\sigma^2 \approx 583$, while for the average $\bar X$, the corresponding parameters are $\mu = 7$ and $\sigma^2 \approx 0.0583$.
One useful feature of this normal approximation is that, for any given deviation of the observed mean from the expected value, we can divide it by the standard deviation $\sigma$ of the distribution (which is simply the square root of the variance $\sigma^2$) and compare this with a precomputed table of confidence intervals, or a simple diagram like the image below, to see how likely the observed mean is to fall that far from the expectation.
Image from Wikimedia Commons by M. W. Toews, based (in concept) on figure by Jeremy Kemp; used under the CC-By 2.5 license.
For example, the average $\bar X$ of 100 rolls of two six-sided dice has a standard deviation of $\sigma = \sqrt{35 \mathbin/ 600} \approx 0.2415$. Thus, an observed mean of 6 or less would be $1/\sigma \approx 4.14 \ge 4$ standard deviations below the expected value of 7, and thus (based on the table on the Wikipedia page I linked to above) occurs with probability less than 0.00006334. A more accurate calculation (using the formula $\frac12\left(1+\operatorname{erf}\left(\frac{x-\mu}{\sigma\sqrt2}\right)\right)$ also given on the Wikipedia page) yields an approximate probability of 0.0000173 of the observed mean falling this far from the expectation, which is fairly close to the value calculated numerically above without using the normal approximation.