I'm going to put in some explanations of things that it sounds like you already know, for the sake of any other readers.

One way of calculating observables in quantum field theory is by performing a path integral. That is, an integral over all configurations of the fields in your system. We start with a partition function

    Z = ∫ D[fields] e^{i S[fields]}

where S is the action, a spacetime integral of the lagrangian,

    S[fields] = ∫ dSpacetime L[fields[spacetime]]

which measures the difference between the kinetic and potential energies.

We can calculate expecation values of an observable O, which is a functional of the fields, by taking averages.

    < O[fields] > = 1/Z ∫ D[fields] e^{i S[fields]} O[fields] 

This integral cannot be done on the cheap with monte carlo, because it's got e^{iS} in it, and you'd need to average a bunch of phases to get a real number. This is a sign problem. Instead, we Wick rotate to Euclidean time and wind up with

    < O[fields] > = 1/Z ∫ D[fields] e^{-S[fields]} O[fields]

This looks much better, and in fact looks like statistical mechanics! Since it is positive-definite, we can treat e^{-S} as a probability distribution. But we still have a problem implementing this computationally, because spacetime is continuous but we have a finite amount of RAM. So we discretize spacetime (typically on a four-dimensional hypercubic lattice) and impose (typically periodic) boundary conditions. This turns the above horrible infinite-dimensional integral into a still-horrible but now only gargantuan-dimensional integral. We replace the integral of S with a discrete sum:

    S_lattice[fields] = ∑[Spacetime] L[fields[spacetime]]

and generate a Markov Chain of field configurations using the Metropolis algorithm for importance sampling. This means we take our current configuration, conf_i and generate a new proposal conf_j. We transition our Markov Chain to conf_j if the weight e^{-S_lattice[conf_j]} is bigger than the weight with conf_i. If the weight is smaller, we transition to conf_j with probability weight_j/weight_i. Otherwise we reject the proposal and conf_i repeats itself in our Markov chain. This algorithm is known to be exact in the large statistics limit. By doing this our Markov chain spends the most time near field configurations with small action that make a large contribution to our integral. Stepping through the Markov chain is very expensive, and every configuration produced is valuable.

This produces an ensemble of N configurations. Because our configurations are generated according to the weight e^{-S} we can estimate simply:

    < O[fields] > = 1/N ∑[field configuration] O[fields] 
                    + error of order 1/sqrt(N)

So, if you let me burn enough electricity I can crank on N and make my uncertainty for the expectation values smaller and smaller. (If you generated according to a flat distribution rather than via importance sampling, you'd need to stick the e^{-S[fields]} in this average.)

Once we have an ensemble of field configurations, we can measure all sorts of observables on them. Typically one measures order parameters or correlation functions, but it depends on what physics you're interested in.

What approximations have I made? First, I calculate on a finite, discrete spacetime volume. Second, I don't have an infinite budget, so I have statistical uncertainty in my result, much like an experimentalist. In fact, we usually say we "measure" the observable O. Of course we're not interrogating nature, but a system of axioms (or something). For our purposes, a computer is like a telescope, but it lets you see into the Platonic world.

To control the continuum limit, you need to redo the above Monte Carlo procedure a few times with different lattice spacings (holding the hadronic spectrum fixed). To control the infinite-volume extrapolation you need to redo the above MC procedure a few times with different numbers of lattice sites but otherwise fixed parameters.

I have neglected some details. For example, there are different ways to discretize, and some ways are better than others (but typically that corresponds to a more expensive method). Depending on your discretization you may have different input parameters you can fiddle with. In QCD, the QCD coupling constant and the quark masses are the ones you can fiddle with. Take some set of parameters and compute, for example, the hadronic spectrum. If the ratios of the masses match those we see in nature, we rejoice. We can then compare the mass computed from the lattice (which is a dimensionless number) and compare it to the masses in the real world (which are measured in GeV). The conversion factor is the lattice spacing (usually on the scale of .1 to 0.05 fm, which we convert to GeV using hbar=1), m_lattice = m_physical * a_lattice. If you are near the continuum limit, anything dimensionful you need to measure on the lattice will simply be converted by the correct powers of the lattice spacing, and any errors will scale like some higher power of the lattice spacing. That this is true is essentially a reflection of the fact that QCD has a negative beta function and becomes weakly coupled / perturbative at large energies / short distances. In theories without a UV fixed point you cannot take a continuum limit---in perturbative language, this is the Landau pole coming to bite you.

As an aside, from this perspective, Wilson's renormalization program is much less mysterious: if you want to hold physics fixed and change your lattice spacing, you'd better change the bare, dimensionless parameters you stuck into your action.