Probability that number of heads exceeds sum of die rolls Let $X$ denote the sum of dots we see in $100$ die rolls, and let $Y$ denote the number of heads in $600$ coin flips. How can I compute $P(X > Y)?$

Intuitively, I don't think there's a nice way to compute the probability; however, I think that we can say $P(X > Y) \approx 1$ since $E(X) = 350$, $E(Y) = 300$, $\text{Var}(X) \approx 292$, $\text{Var}(Y) = 150$, which means that the standard deviations are pretty small.
Is there a better way to approach this problem? My explanation seems pretty hand-wavy, and I'd like to understand a better approach.
 A: The exact answer is easy enough to compute numerically — no simulation needed.  For educational purposes, here's an elementary Python 3 script to do so, using no premade statistical libraries.
from collections import defaultdict

# define the distributions of a single coin and die
coin = tuple((i, 1/2) for i in (0, 1))
die = tuple((i, 1/6) for i in (1, 2, 3, 4, 5, 6))

# a simple function to compute the sum of two random variables
def add_rv(a, b):
  sum = defaultdict(float)
  for i, p in a:
    for j, q in b:
      sum[i + j] += p * q
  return tuple(sum.items())

# compute the sums of 600 coins and 100 dice
coin_sum = dice_sum = ((0, 1),)
for _ in range(600): coin_sum = add_rv(coin_sum, coin)
for _ in range(100): dice_sum = add_rv(dice_sum, die)

# calculate the probability of the dice sum being higher
prob = 0
for i, p in dice_sum:
  for j, q in coin_sum:
    if i > j: prob += p * q

print("probability of 100 dice summing to more than 600 coins = %.10f" % prob)

Try it online!
The script above represents a discrete probability distribution as a list of (value, probability) pairs, and uses a simple pair of nested loops to compute the distribution of the sum of two random variables (iterating over all possible values of each of the summands).  This is not necessarily the most efficient possible representation, but it's easy to work with and more than fast enough for this purpose.
(FWIW, this representation of probability distributions is also compatible with the collection of utility functions for modelling more complex dice rolls that I wrote for a post on our sister site a while ago.)

Of course, there are also domain-specific libraries and even entire programming languages for calculations like this.  Using one such online tool, called AnyDice, the same calculation can be written much more compactly:
X: 100d6
Y: 600d{0,1}
output X > Y named "1 if X > Y, else 0"

Under the hood, I believe AnyDice calculates the result pretty much like my Python script does, except maybe with slightly more optimizations.  In any case, both give the same probability of 0.9907902497 for the sum of the dice being greater than the number of heads.
If you want, AnyDice can also plot the distributions of the two sums for you.  To get similar plots out of the Python code, you'd have to feed the dice_sum and coin_sum lists into a graph plotting library like pyplot.
A: It is possible to do exact calculations.  For example in R
rolls <- 100
flips <- 600
ddice <- rep(1/6, 6)
for (n in 2:rolls){
  ddice <- (c(0,ddice,0,0,0,0,0)+c(0,0,ddice,0,0,0,0)+c(0,0,0,ddice,0,0,0)+
            c(0,0,0,0,ddice,0,0)+c(0,0,0,0,0,ddice,0)+c(0,0,0,0,0,0,ddice))/6}
sum(ddice * (1-pbinom(1:flips, flips, 1/2))) # probability coins more
# 0.00809003
sum(ddice * dbinom(1:flips, flips, 1/2))     # probability equality
# 0.00111972
sum(ddice * pbinom(0:(flips-1), flips, 1/2)) # probability dice more
# 0.99079025

with this last figure matching BruceET's simulation
The interesting parts of the probability mass functions look like this (coin flips in red, dice totals in blue)

A: Another way is by simulating a million match-offs between $X$ and $Y$
to approximate $P(X > Y) = 0.9907\pm 0.0002.$ [Simulation in R.]
set.seed(825)
d = replicate(10^6, sum(sample(1:6,100,rep=T))- 
       rbinom(1,600,.5))
mean(d > 0)
[1] 0.990736
2*sd(d > 0)/1000
[1] 0.0001916057   # aprx 95% margin of simulation error


Notes per @AntoniParellada's Comment:
In R, the function sample(1:6, 100, rep=T) simulates 100 rolls a fair die;
the sum of this simulates $X$. Also rbinom is R code for simulating
a binomial random variable; here it's $Y.$ The difference is $D = X - Y.$
The procedure replicate makes a vector of a million differences d.
Then (d > 0) is a logical vector of a million TRUEs and FALSEs, the mean of which is its proportion of TRUEs--our Answer. Finally, the last statement
gives the margin of error of a 95% confidence interval of the proportion
of TRUEs (using 2 instead of 1.96), as a reality check on the accuracy
of the simulated Answer. [With a million iterations one ordinarily expects
2 or 3 decimal paces of accuracy for probabilities--sometimes more for
probabilities so far from 1/2.]
A: A bit more precise:
The variance of a sum or difference of two independent random variables is the sum of their variances. So, you have a distribution with a mean equal to $50$ and standard deviation $\sqrt{292 + 150} \approx 21$.
If we want to know how often we expect this variable to be below 0, we can try to approximate our difference by a normal distribution, and we need to look up the $z$-score for $z = \frac{50}{21} \approx 2.38$.
Of course, our actual distribution will be a bit wider (since we convolve a binomial pdf with a uniform distribution pdf), but hopefully this will not be too inaccurate.
The probability that our total will be positive, according to a $z$-score table, is about $0.992$.
I ran a quick experiment in Python, running 10000 iterations, and I got $\frac{9923}{10000}$ positives. Not too far off.
My code:
import numpy as np
c = np.random.randint(0, 2, size = (10000, 100, 6)).sum(axis=-1)
d = np.random.randint(1, 7, size = (10000, 100))
(d.sum(axis=-1) > c.sum(axis=-1)).sum()
--> 9923

