Monte Carlo modeling takes random samples of input distributions. You put  together a behavioral model (for example, a thermal network, or a linearization of data from a physical or numerical experiment). The output of the behavioral model is a distribution of the output parameter. (Thanks to Garron Morris and CoolingZone for the helpful articles.)

As an example (not the one I did for a client...), consider a "simple" lidded flip chip package with a heat sink on it. The average impinging velocity has some variation from part to part  within the system, and possibly also system to system. The performance of the two thermal interface materials also varies.

Using a spreadsheet, I extracted a fast behavioral model for the thermal resistance of the heat sink in impinging flow (based on a paper). The maximum die temperature is a function of this thermal resistance, and also of the power (which varies part to part even under constant load).

The input distributions, which in this case are fictitious (note zero mean values for TIM distributions -- I rolled nominal values into the theta equation), might look something like this:

Using these input distributions, you get the thermal resistance and chip temperature distributions looking like this -- or at least, these are examples of how the distributions might look (you're sampling, remember -- it takes a lot of samples to get repeatable distributions).

The solution that looks like it'll work fine, based on nominal values, to keep the maximum chip temperature below 100 C will  actually produce quite a few higher chip temperatures at the specified ambient temperature (which you could also use as an input random variable, by the way -- if you could get better information about your  system's install base).