Puzzled about the non-teaching of Monte Carlo method(s) for error analysis

[FactChecker]

Since R is always positive, it was a given that the distribution is not Gaussian. There are some common distributions that are always positive. The log-normal distribution might be a good fit with the correct parameters. It would be interesting to calculate the mean and variance of the associated ln(R) values and see if the log-normal distribution with those parameters is a good fit. (see https://en.wikipedia.org/wiki/Log-normal_distribution#Generation_and_parameters)

 

[hutchphd]

It turns out I didn't make any mistake there (I think!). It may seem counter intuitive, so here is the same plot (well another such simulation) with the mean, 68.27 interval, 95.45 and 99.73 intervals (1, 2 and 3 SD for a Gaussian). It may be an artifact of the histogram, which contains "only" 250 bins, so beware.

There is no way 68% of the curve area is between the orange lines.
What is a count??
Of course I have no way to check this as previously discussed.
The one SD number will be ~10% of the average=30ohms. You have done something wrong.
For a Guassian FWHM (Full width half max) is 2.3 sigma

 

[hutchphd]

Since R is always positive, it was a given that the distribution is not Gaussian.

Of course that is true. But for things within one SD it will be a very good approximation. The "wings" will be less so. But unless I misunderstand this calculation is nonsense.

 

[Vanadium 50]

There is no way 68% of the curve area is between the orange lines.

I agree. The first blue plot has a FWHM of around 75-80Ω, which means a standard deviation around 32-34Ω, which is 10-11%. That agrees with what was put in.

+16 -14 ohm cannot possibly be right. If your worst measurement is 10%, you can't get a 5% measurement by combining it with other uncertainties. It doesn't match the plot.

Somehow this is evidence that the analytic formula is wrong? I would say instead that it's evidence that an over-reliance on MC makes it harder to spot even obvious errors.

 

[FactChecker]

I would say instead that it's evidence that an over-reliance on MC makes it harder to spot even obvious errors.

Overconfidence is probably a better word than over-reliance (although the OP is actually questioning the results). Results should be critically examined, whether they are from classical analysis or from MC simulations. In fact, this resistance example is as simple as you will ever see. (I think that the problem statement that current and voltage are assumed to be uncorrelated might invalidate any approach.) When the subjects become network theory, resource allocation, spares positioning, queueing theory, random walk, etc., even simple problems are difficult analytically. In those, the slightest complication can make analytical approaches unreasonable.
I am not saying that one should not thoroughly understand the statistical analysis methods. I am just saying that Monte Carlo simulations are the only practical way to address a great many problems. I believe that the OP implied that the analytical theory should be de-emphasized. I disagree with that. I also disagree with comments implying that Monte Carlo methods are just an easy and sloppy substitute for good analysis. MC methods are an essential addition to the toolset.

 

[Vanadium 50]

Overconfidence is probably a better word than over-reliance

That's fair.

I am not saying that one should not thoroughly understand the statistical analysis methods. I am just saying that Monte Carlo simulations are the only practical way to address a great many problems

I agree with that (example to follow). I don't think anyone is arguing against that.

I believe that the OP implied that the analytical theory should be de-emphasized.

That's what people are arguing against.

I have no objection to MC methods. Use 'em myself. But it's a case of the right tool for the right job.

Promised example: I have a copper bus bar of dimensions (l,w,h) and with (possibly different) uncertainties on each dimension. However, the total uncertainty on the volume is much less than any of those uncertainties. (Perhaps it was very accurately weighed). What is the uncertainty on its resistance?

While I think I could solve it analytically (but only because the cross-sectional shape factors out), it's probably simpler to run an MC.

 

[fluidistic]

I am sorry, I am destroyedly tired. I also want to understand what's wrong with the data or histogram.
Python code that generates a plot like I've provided:

import random as rdn
import statistics as sts
import matplotlib.pyplot as plt
import numpy as np

def one_run():
    I = rdn.gauss(0.025e-3, 0.1*0.025e-3)
    V = rdn.gauss(77e-4, 0.02*77e-4)
    R = V/I
    return R

resistances = []
for run in range(100000):
    result = one_run()
    resistances.append(result)

the_median = sts.median(resistances)
the_mean = sts.mean(resistances)
lower_confidence = np.percentile(resistances, 100-68.3)
upper_confidence = np.percentile(resistances, 68.3)
lower_2SD = np.percentile(resistances, 4.65)
lower_3SD = np.percentile(resistances, 0.27)
upper_2SD = np.percentile(resistances, 95.45)
upper_3SD = np.percentile(resistances, 99.73)

print('median', the_median)
print('mean', the_mean)
print('minus window', the_median - lower_confidence)
print('upper window', upper_confidence - the_median)
#print(resistances)
fig, ax = plt.subplots(1, 1)

ax.hist(resistances, bins=250)
plt.axvline(lower_2SD, color='yellow')
plt.axvline(upper_2SD, color='yellow')
plt.axvline(lower_3SD, color='green')
plt.axvline(upper_3SD, color='green')
plt.axvline(the_median, color='red')
plt.axvline(lower_confidence, color='orange')
plt.axvline(upper_confidence, color='orange')

ax.set_xlabel('Resistance (ohm)')
ax.set_ylabel('Count')
#plt.show()
plt.savefig('R_hist.jpg')
plt.close()

Of course, I do not ask for a debug, but if you want to check it out yourself, feel free to use the code. I will debug when I can, as I would like to understand what is going on.

Edit: I figured out the mistake, sorry. (Facepalm, shameful face). I will correct it. The histogram is correct I think. But not the error bars, as you have pointed out.

 

[fluidistic]

I am sorry, I am destroyedly tired. I also want to understand what's wrong with the data or histogram.
Python code that generates a plot like I've provided:

import random as rdn
import statistics as sts
import matplotlib.pyplot as plt
import numpy as np

def one_run():
    I = rdn.gauss(0.025e-3, 0.1*0.025e-3)
    V = rdn.gauss(77e-4, 0.02*77e-4)
    R = V/I
    return R

resistances = []
for run in range(100000):
    result = one_run()
    resistances.append(result)

the_median = sts.median(resistances)
the_mean = sts.mean(resistances)
lower_confidence = np.percentile(resistances, 100-68.3)
upper_confidence = np.percentile(resistances, 68.3)
lower_2SD = np.percentile(resistances, 4.65)
lower_3SD = np.percentile(resistances, 0.27)
upper_2SD = np.percentile(resistances, 95.45)
upper_3SD = np.percentile(resistances, 99.73)

print('median', the_median)
print('mean', the_mean)
print('minus window', the_median - lower_confidence)
print('upper window', upper_confidence - the_median)
#print(resistances)
fig, ax = plt.subplots(1, 1)

ax.hist(resistances, bins=250)
plt.axvline(lower_2SD, color='yellow')
plt.axvline(upper_2SD, color='yellow')
plt.axvline(lower_3SD, color='green')
plt.axvline(upper_3SD, color='green')
plt.axvline(the_median, color='red')
plt.axvline(lower_confidence, color='orange')
plt.axvline(upper_confidence, color='orange')

ax.set_xlabel('Resistance (ohm)')
ax.set_ylabel('Count')
#plt.show()
plt.savefig('R_hist.jpg')
plt.close()

Of course, I do not ask for a debug, but if you want to check it out yourself, feel free to use the code. I will debug when I can, as I would like to understand what is going on.

Edit: I figured out the mistake, sorry. (Facepalm, shameful face). I will correct it. The histogram is correct I think. But not the error bars, as you have pointed out.

Tadum:
Uncertainties that doesn't vary much from the Gaussian (which makes sense, assuming the histogram has no hidden surprises).
minus window 28.76204359267797
upper window 34.91842045829611
Updated plot

 

[hutchphd]

This seems a perfect cautionary tale to me. .

 

[Vanadium 50]

minus window 28.76204359267797
upper window 34.91842045829611

If you can't be accurate, at least be precise...

 

[FactChecker]

This seems a perfect cautionary tale to me.

It looks like the Monte Carlo simulation was only a couple of simple lines of code that were always correct. It was the data post-analysis that went wrong. The real cautionary tale might about a broader issue: The assumed zero correlation of the voltage and current through a (fixed?) resistance. That is something to be concerned about no matter what approach is taken. IMHO, there can be less emphasis on difficult theoretical analysis and more on frequent, fundamental, issues like whether variables are correlated, self-selecting samples, and other sources of bias that are so often overlooked.

 

[fluidistic]

It looks like the Monte Carlo simulation was only a couple of simple lines of code that were always correct. It was the data post-analysis that went wrong. The real cautionary tale might about a broader issue: The assumed zero correlation of the voltage and current through a (fixed?) resistance. That is something to be concerned about no matter what approach is taken. IMHO, there can be less emphasis on difficult theoretical analysis and more on frequent, fundamental, issues like whether variables are correlated, self-selecting samples, and other sources of bias that are so often overlooked.

Sure, one must be cautious. Also, although I have performed a simple sensitivity test, which showed the effect of improving the accuracy of either I or V and concluded like the analytical formula, I haven't obtained numbers but just histograms. I should investigate how I could objectively determine the improvements.
About the non correlation between I and V, do you mean that in the real world, R is a function of temperature and that using a current creates a Joule effect that affects T, thus R and thus V? Or do you have something else in mind?

 

[FactChecker]

About the non correlation between I and V, do you mean that in the real world, R is a function of temperature and that using a current creates a Joule effect that affects T, thus R and thus V? Or do you have something else in mind?

Sorry. I may have missed some explanation. It is not clear to me what the experiment is. If a power supply is placing a voltage across a single resistor with constant resistance, the voltage and current are definitely correlated. You are simulating uncorrelated voltage and current. Are you talking about measurements with uncorrelated measurement errors? Are you talking about real voltage errors in the power supply? Both? Something else? I am confused.

 

[fluidistic]

Sorry. I may have missed some explanation. It is not clear to me what the experiment is. If a power supply is placing a voltage across a single resistor with constant resistance, the voltage and current are definitely correlated. You are simulating uncorrelated voltage and current. Are you talking about measurements with uncorrelated measurement errors? Are you talking about real voltage errors in the power supply? Both? Something else? I am confused.

There's no real experiment, so it's up to us. I mentioned a current source (not a voltage source), but we're not forced to focus on this. I am just wondering the sources of correlations between I and V in a real life experiment.
For the simulations I attempted to follow hutchphd's example where instead of using an ammeter with 10 per cent error, I assumed that the current source has a 10 percent error on what the current is set at. Then a voltmeter with 2 percent accuracy is measuring the potential drop across the resistor.

 

[FactChecker]

There's no real experiment, so it's up to us. I mentioned a current source (not a voltage source), but we're not forced to focus on this. I am just wondering the sources of correlations between I and V in a real life experiment.
For the simulations I attempted to follow hutchphd's example where instead of using an ammeter with 10 per cent error, I assumed that the current source has a 10 percent error on what the current is set at. Then a voltmeter with 2 percent accuracy is measuring the potential drop across the resistor.

Thanks. That answers my question. So you are measuring a fixed, but unknown resistance using meters that have the errors you describe. The voltage and current errors are simply measurement errors that are independent of each other. That makes sense to me now. So your analysis and Monti Carlo simulation are about the distribution of the calculated resistance using Ohm's law with the measurement errors. I guess that the means of the MC voltage and current simulations are set for the specific true resistance. I would have to question whether the mean of a ratio is simply the ratio of the means.

 

[Stephen Tashi]

I would have to question whether the mean of a ratio is simply the ratio of the means.

Also, the ratio of two independent normally distributed random variables with zero means has a Cauchy distribution, whose mean does not exist. (I gather the idea is to estimate the true resistance as the mean of the V/I values.) The ratio of independent normals with non-zero means may be better behaved.

A paper by Marsaglia about dealing with the ratio of normal random variables is:

https://www.google.com/url?sa=t&rct...sg=AFQjCNEgO1dvktreWiL-rt-ZPcS3K1FmYQ&cad=rja

As to the merits of pencil-and-paper methods of "error propagation", my impression (based only on questions that appear in the forum about error propagation) is that (statistically!) instruction in the pencil and paper methods seems to teach a set of procedures, but does not connect those procedures with the theory of probability models. Many questioners ask about how to compute "uncertainties" in things and are confident that "uncertainty" has an obvious meaning. It may indeed have a well defined meaning in their particular technical discipline, but I don't detect this meaning is identical to a particular concept in probability and statistics. Of course, there is a potential bias in my sample since people with less understanding of the subject would be more inclined to ask questions about it.

 

[fluidistic]

Thanks. That answers my question. So you are measuring a fixed, but unknown resistance using meters that have the errors you describe. The voltage and current errors are simply measurement errors that are independent of each other. That makes sense to me now. So your analysis and Monti Carlo simulation are about the distribution of the calculated resistance using Ohm's law with the measurement errors. I guess that the means of the MC voltage and current simulations are set for the specific true resistance. I would have to question whether the mean of a ratio is simply the ratio of the means.

I think the median nails it pretty accurately. The mean is going to be biased towards greater values. Part of the reason why I asked why the GUM states to report the mean and SD rather than the median and confidence intervals.

 

[Stephen Tashi]

Part of the reason why I asked why the GUM states to report the mean and SD rather than the median and confidence intervals.

A simple argument for reporting the mean and SD:
Suppose we do N independent replications of a simulation and each replication produces one value of Y. Then, by the central limit theorem, the mean value of all the Y's is approximately normally distributed for large N , provided the distribution of Y (as an individual random variable) has finite variance. (The distribution of Y need not be a normal distribution for this to hold.) So if the purpose of doing the simulation is to estimate the "true" value Y by using the mean of the Y data, the estimated standard deviation of Y is relevant to computing confidence intervals.

In the above approach, the size of confidence intervals depends on the number of samples (replications) since that number enters into the calculation of the sample standard deviation of the Y data. By contrast, you suggest looking at a histogram and its median to compute a confidence interval. How does the size of this confidence interval depend on the number of replications you did to create the histogram? If the number you report is not a function of the number of replications of the simulation, then can it be interpreted as a "confidence interval" or should we use some other terminology for it?

 

[hutchphd]

Also, the ratio of two independent normally distributed random variables with zero means has a Cauchy distribution, whose mean does not exist.

These variables do not have zero means

 

[FactChecker]

The distribution of the ratio of uncorrelated normal variables with non-zero means is discussed here. They give the equation for the variance. It uses the individual means and variances. I do not see anything about the mean. It makes me tired.
This is something where I would feel more comfortable with the Monte Carlo simulation than with my analysis, and it really is a very basic thing.

 

[hutchphd]

Zero over zero is never fun.

 

[Stephen Tashi]

This is something where I would feel more comfortable with the Monte Carlo simulation than with my analysis, and it really is a very basic thing.

And the most basic thing, is: What is the problem being analyzed? The question highlights the distinction between the usual type of "confidence interval" analysis versus "error propagation" analysis.

In confidence interval analysis, we usually have data in the form of independent samples. The confidence interval analysis depends on the data and the amount of data we have.

By contrast, the example being simulated assumes some given distributions. There is no mention of how the parameters of these distributions relate to any data. It is as if the parameters of the distributions are handed to us and assumed to be exact.

The general form of confidence interval analysis is that we have a formula $\hat{c}$ for estimating some parameter $c$ of a distribution as a function of data. The usual form of a confidence interval is $( \hat{c} - L/2, \hat{c} + L/2)$ where $L$ is positive constant. A confidence interval depends on data that has randomness, so it is a randomly generated interval. The probability associated with confidence intervals of length $L$ is the probability that the true value of the parameter lies within at least one of these randomly generated intervals. (This probability is not associated with single confidence intervals. For example, a "95 % confidence" is not a probability that can be associated with a single interval like $(25.3 - L/2, 25.3 + L/2)$).

So how do the simulation results relate to a confidence interval analysis? I think the histogram of $V/I$ is relevant to analyzing whether a single measurement (a sample of size 1) of $(V,I)$ is within a certain distance of the true value of $R$. If multiple measurements are taken, further analysis is needed. However, the parameters of the distributions for $V,I$ used in the simulation results are not presented as estimates from data. It appears logically contradictory to analyze the situation of multiple $V/I$ measurements using a simulation that ignores them.

 

[FactChecker]

It's a very simple Monte Carlo simulation where the V and I measurements are modeled as Gaussian. I think that the sample results can be used for probability distributions and confidence intervals as with any other random sample. The validity of the model to represent the particular resistance problem might be questionable.

 

[Stephen Tashi]

It's a very simple Monte Carlo simulation where the V and I measurements are modeled as Gaussian. I think that the sample results can be used for probability distributions and confidence intervals as with any other random sample.

I agree that the simulation is simple. But what is the analysis? For example, what is the length of a 95% confidence interval centered on the mean ( or median if we prefer) of 10 pairs of $V/I$ measurements?

 

[FactChecker]

I guess I have to take back any strong statement about the confidence intervals. I would treat the Monte Carlo sample the same as any other statistical sample.
I am also troubled by an analytical approach that allows it to be a ratio of normal random variables. Suppose we ignore the problem of division by zero. What do we say if the voltage measurement is positive and the current measurement is negative? Or vice versa. Do we just throw out that possibility? And what if both happen to be negative? Since they are assumed to be uncorrelated, that would be possible. Would we treat that as more valid just because they are not mixed positive/negative and they give a positive resistance? This resistor problem is a hypothetical problem and I am not sure that we are modeling it properly.
The problem hurts my head. I think I will leave this problem to others.

 

[Vanadium 50]

has a Cauchy distribution

You mean a Breit-Wigner, right? Or maybe a Lorentzian?

 

[hutchphd]

The problem hurts my head. I think I will leave this problem to others.

I am certainly glad I was never "educated" in this stuff. And I used it a lot with great practical success, largely oblivious I guess.
All I require is that the deviations (rms errors) be small compared means, and they be "uncorrelated". Then the leading order terms for the rms deviations are as described independent of the form of the distributions. Of course if the Taylor expansion for the functional dependence blows up there is a problem, but in the real world this seldom happens.
These techniques are extraordinarily useful and robust. In my experience the only places requiring some care are low probability events (the wings of the distribution) where wrong assumptions will bite you. Do not be afraid. Say the magic words: "central limit theorem".

 

[Vanadium 50]

First, the thing we really want is the probability that the true value is inside the error bars to be 68%. That's not well-defined. (And to also require that the probability that the true value is also inside twice the error bars to be 90% is even less well-defined). So our whole statistical technique, analytic or Monte Carlo, isn't built on a mathematically rigorous foundation. But it's the best we have. And "not mathematically rigorous" is not the same as useless - there's real meaning in comparing error bars of 1%, 5% and 10%, even if none of these are exactly what we hoped to know.

What do we say if the voltage measurement is positive and the current measurement is negative?

Here's where you have to think like a physicist, not a mathematician. If my knowledge of the current is so poor I can't tell which way it is flowing, or even if it is flowing at all, I shouldn't be using it to calculate the resistance.

 

[FactChecker]

In my experience the only places requiring some care are low probability events (the wings of the distribution) where wrong assumptions will bite you. Do not be afraid. Say the magic words: "central limit theorem".

You are correct, of course. I was thinking of the 10% voltage error as being a huge error without realizing that a negative voltage would be 10 standard deviations below the mean. Stranger things have happened, but not since Moses.

 

[Stephen Tashi]

If the goal of the simulation example was to compare the simulation technique with a theoretical pencil-and-paper technique then we'd need to see the pencil-and-paper technique worked out to make the comparison. But to see the pencil-and-paper technique worked out, we'd need to define what problem is being analyzed in the first place!

 

[Dr Transport]

The bottom line is you need to know what the hell you are doing...The Monte Carlo methods give you numbers and very little insight. If you don't understand the theory, you will not know what you are doing.

My PhD advisor said exactly the same thing to me 25 years ago.

 

[FactChecker]

Insight and theory are great when they are correct. But things often get very complicated and confusing. If nothing else, Monte Carlo simulations can often be used to verify theoretical results or to point out errors.

Here is an example showing the difficulty of analyzing a fairly simple queueing problem: https://www.physicsforums.com/threads/waiting-time-in-a-queue-using-poisson-arrival.902175
And the only way I would feel confident of the analytical results is if there was a simulation that supported them. The only way to do real work in queueing theory is with MC simulation. IMHO, the only queueing problems that can be solved analytically are trivial ones.

Here is a problem involving a dice game where the analytical solution is messy and the MC simulation is simple. https://www.physicsforums.com/threads/probability-in-a-dice-game.989492/
Again, I would not be confident of the analytical solution at all if there was not an MC simulation result to support it.

 

[mpresic3]

I had a professor who taught Classical Electrodynamicsi at the graduate level, not from Jackson,but with his own notes involving exclusively, differential forms. I knew a professor that was proposing to teach graduate Classical Mechanics, not from Goldstein, but using category theory, alone. Both suggested these areas were sadly lacking in graduate physics education. After attending the Electrodynamics course, I had (two) graduate courses in Jackson. (I went to graduate school twice).

It may seem I am favor of conventional treatment of the physics curriculum. I think there is danger in uniformity, and it is good for some students to have different tools in their toolbox. However, it is hard to come up with areas in the tight physics curriculum that could be left out. Certainly, including MC methods at the expense of other important topics is going to be objected too by others.

It seems like, when the poster becomes the instructor head of the course, he or she can then teach whatever he or she wants. There is quite a bit academic freedom in the USA anyway. The professor who taught differential forms did not get much pushback. The professor who proposed category theory (as far as I know) did not get his course, because no student was interested in taking the course.

Also, the training of a physicist contains more than just physics courses. Physicists can run into MC techniques in computer science courses, or statistics courses. A good argument could be made that statistics and probability should be required. Maybe some would say, substitute probability, for Complex Analysis. However, a look at most graduate physics program, seems to regard complex analysis over probability. As I wrote, you can always find somebody to find something they feel should be part of the education, that is overlooked.