Chronology Current Month Current Thread Current Date
[Year List] [Month List (current year)] [Date Index] [Thread Index] [Thread Prev] [Thread Next] [Date Prev] [Date Next]

[Phys-l] monte carlo uncertainty analysis



Hello everyone,

I was reading John Denker's excellent document on uncertainties, http://www.av8n.com/physics/uncertainty.htm, and thought it'd be fun to program the Monte Carlo simulation technique he describes in it. Although a spreadsheet can be used, and is given on that site, I find it easier to use something like Python for analysis. In my code, you specify the value and the uncertainty, like:

b=Quantity(0.4,0.04)

or

d=Quantity('0.4 +- 10%')

and some other ways.

You then then ask it to evaluate some expression, like:

x=evaluate('(((a+b)+c)+d)+e')

It then samples (default 10000 samples) all of the values you specified, with a Gaussian distribution with mean and standard deviation set by the value and uncertainty of the quantity. It returns then the mean and standard deviation of the result.

I use as demonstrations two examples from John Denker's site. The following code:

a=Quantity(2) # no uncertainty
b=Quantity(0.4,0.04) # 10 percent
c=Quantity('0.4 +- 0.04') # 10 percent
d=Quantity('0.4 +- 10%') # 10 percent
e=b

x=evaluate('(((a+b)+c)+d)+e')

gives

x = 3.5984327533 +- 0.0804616538054


and the following:

Mg24_mass=Quantity("23.9850423(8)")
Mg25_mass=Quantity("24.9858374(8)")
Mg26_mass=Quantity("25.9825937(8)")

Mg24_abundance=Quantity("78.99(4)")
Mg25_abundance=Quantity("10.01(1)")
Mg26_abundance=Quantity("11.01(3)")

x=evaluate("""Mg26_mass*Mg26_abundance/100+
Mg25_mass*Mg25_abundance/100+
Mg24_mass*(100-Mg26_abundance-Mg25_abundance)/100""")

gives

Mg mass with sum rule: 24.3051502306 +- 0.000603339566833

and (for no sum rule)


x=evaluate("""Mg26_mass*Mg26_abundance/100+
Mg25_mass*Mg25_abundance/100+
Mg24_mass*Mg24_abundance/100""")


Mg mass without sum rule: 24.3076501938 +- 0.0127462356118

You can get a copy of the code from my projects link off of my website, http://web.bryant.edu/~bblais.

I'd love to hear from anyone interested in this, or if there are any additions/suggestions.


thanks,

Brian Blais

--
Brian Blais
bblais@bryant.edu
http://web.bryant.edu/~bblais