What is the fundamental reasoning for anti aliasing using multiple random samples within a pixel? In graphics, it's common to take multiple samples within the bounds of a pixel and combine them together (most commonly just doing an average) for a final sample pixel color.  This has the effect of anti aliasing an image.

On one hand this makes sense to me because what you are effectively doing is integrating the color of the pixel over the area that the pixel represents.  In this line of thinking, averaging "random" samples seems to be the ideal setup, for doing monte carlo integration.  ("random" could be stratified, blue noise based, low discrepency sequences etc)

On the other hand, this feels wrong (or at least not as correct as it could be) from a digital signal processing point of view.  From that point of view, it feels like we are taking a lot of samples and then downsampling using a box filter (box blur) to get the final pixel value.  In that light, it seems like the ideal thing to do would be to use sinc filtering instead of averaging the samples.  I could see that the box filter is a cheaper aproximation to sinc thinking along these lines.

This leaves me a bit confused.  Is the core idea that we are integrating the pixel area and averaging is correct?  Or is it that we are downsampling and should be using sinc, but are using a box filter because it's fast?

Or is it something else entirely?

A little bit related: Anti-aliasing / Filtering in Ray Tracing
From a signal processing point of view, you're sampling a continuous-domain signal, and you need to filter it to get rid of frequencies beyond the Nyquist limit. It's that filtering that leads to integrating over the pixel area—or more generally, integrating over the support of your antialiasing kernel (which need not be a box).

Consider your rendering function that takes a sample point $x, y$ in screen space and returns back the color found at that point. (Let's ignore any issues of random sampling for the moment, and assume it returns a "perfect" rendered color for that specific point.) This function effectively defines a 2D continuous-domain signal. Or to put it another way, it defines an infinite-resolution image, since nothing prevents this function from having features at arbitrarily small scales. In terms of frequency domain: the function is not band-limited; it can include components of arbitrarily high spatial frequencies.

Now you want to convert it down to a finite number of pixels. Much like digitizing an audio signal, when you sample it, you'll get aliasing unless you first eliminate frequencies beyond the Nyquist limit imposed by the sampling rate. In other words, you have to get rid of features smaller than the pixel grid. To do this, you apply a low-pass filter. The ideal low-pass filter is the sinc function, but for various reasons of practicality we use other filters (that don't perfectly eliminate frequencies beyond the Nyquist limit, but they at least attenuate them).

Low-pass filtering is done by convolution. If $f(x, y)$ is the rendering function and $k(x, y)$ is a filter kernel, then we can define the low-pass filtered image as:

$$f_\text{filtered}(x, y) = \iint f(x', y') \, k(x' - x, y'- y) \, dx' \, dy'$$

Then the image is safe to sample, so the final pixel values are gotten by just evaluating $f_\text{filtered}$ at the pixel coordinates.

If $k$ is a box filter, which looks like $k = 1$ within the pixel box and $k = 0$ elsewhere, then this simplifies to just integrating $f$ over the pixel box. But as noted, the box filter isn't so great and there are better choices such as tent, bicubic, and Gaussian filters.

Anyway, now that we have an integral, we can use Monte Carlo for it and combine it with any other integrals we might like to do—for lighting, motion blur, and so on. We can even apply importance sampling to the $k$ factor in the integral, by generating samples for each pixel that are distributed around the pixel center according to $k$.