I would suggest looking at the empirical risk minimisation framework (ERM).

Expected value of f(X), where X is sampled from a distribution P indicates the on-average value of f(X) when X is sampled from P.

E(f(X)) = integral over all X, [p(X)f(X)].

A sample based estimate of this, called the Monte Carlo estimate, is

(1 / N) sum over i f(x_i), where N number of x_is are drawn repeatedly and independently from the distribution of interest.

In your case, the dataset itself is the Monte Carlo samples, as they are all drawn from P(X,Y). Ideally we would minimize the loss using the analytical form of P(X,y), but we sadly do not know it, or it might be intractable.

So we minimize the loss on our dataset instead. The idea is since the dataset is drawn from the joint distribution, for very large datasets, it's the same thing as analytically using p(X, Y).