The Ising model can be solved exactly in one dimension, and in two with no field. In three, nobody has solved it. If it cannot be solved, counting might do instead. It cannot.
There are 2^N arrangements of N spins. The figure is how that grows: a mere 10 \times 10 lattice already gives 2^{100}, more than the atoms in the universe, and no rate of counting will finish it.
But what we want is an average, and an average does not require counting everything. It is enough to pick out some arrangements, starting with the ones that happen most. The problem is how to pick those out.
We make a rule for moving from arrangement to arrangement. Run it long enough and we want the share of time spent in arrangement a to be proportional to e^{-\beta E_a}. The rule is chosen so that it is.
A sufficient condition for that is P(a)W(a \to b) = P(b)W(b \to a): the flow from a to b balances the flow back from b to a, for every pair. This is called detailed balance.
The two arrows in the figure are those flows, drawn with thickness for size. A thin arrow leaves the arrangement with the larger share, a thick one the smaller. Once they balance, the shares stop moving. That is the distribution we were after.
Many rules satisfy detailed balance; the simplest is the Metropolis acceptance. Pick one spin, turn it over, and look at the change in energy \Delta E.
If it falls, always accept. If it rises, accept with probability e^{-\beta\Delta E}. Written together, \min(1, e^{-\beta\Delta E}). The figure is that shape.
When it is cold, moves upward are almost never accepted; when it is hot, almost always. Temperature enters at this one place only, and everything else follows from it. Accepting only the moves that lower the energy would stop the system in its lowest arrangement, which is not the average at any finite temperature. Accepting upward moves is what makes it one.
The figure starts from a scattered arrangement and swaps spins by that acceptance, one after another, at a temperature below the critical one.
At first up and down are finely mixed. As it runs, regions of one direction grow, their boundaries move, and one side eventually spreads. What mean field gave us as falling one way or the other happens here in plain sight.
Taken one move at a time, moves that raise the energy are being accepted too. And still the whole falls. It is the weight e^{-\beta\Delta E} that makes the difference between the two directions.
The figure is the energy over that run. It falls sharply at first, and past some point it stops falling. Where it stops falling is equilibrium.
Even in equilibrium the line does not go still. It keeps moving finely up and down. That is the fluctuation, and its size ties to the heat capacity, exactly as in the piece on how much the energy fluctuates. We can see it here because this method runs a system of finite size.
The part before it settles is not equilibrium, so it is thrown away before averaging. How much to throw away is something the method will not tell you itself. Even so, the answers to models nobody can solve are now available as numbers. The \beta = 0.33 of the three-dimensional Ising model was obtained this way.