Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Mathematica Numerical vs. Monte-Carlo Simulations

  1. May 2, 2016 #1
    I have an integration that doesn't have a solution in the table of integrals. So, I evaluated it using Mathematica using the command NIntegrate. However, when I compare the result with Monte-Carlo simulations, there is a very small constant gap between the two curves. Is it because of the numerical integration accuracy?
     
  2. jcsd
  3. May 2, 2016 #2
    How can I upload an image from my computer?
     
  4. May 2, 2016 #3

    Dale

    Staff: Mentor

    Usually NIntegrate is pretty reliable. The algorithm selection mechanism is very good at choosing the right algorithm for the particular integrand. I would be much more skeptical of the Monte Carlo algorithm. I would look for errors there first.
     
  5. May 2, 2016 #4
    Error like what? In Monte-Carlo simulations, I basically generate random variables, insert them in an expression, and then average the expression over the number of samples to get the average value of the expression.
     
  6. May 2, 2016 #5
    See here, please.
     
  7. May 2, 2016 #6

    Dale

    Staff: Mentor

    Could be anything. You could have a bad random number generator, insufficient number of samples, an error in the arithmetic, numerical precision problems , or any number of other mistakes.

    The point is that it is hard to find a more well tested numerical integrator than NIntegrate. If your problem can be expressed as an integral and evaluated with NIntegrate then that is your answer, the Monte Carlo is just an approximation with a much less sophisticated process for controlling the errors.
     
  8. May 3, 2016 #7

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    Your link goes to an integral from 0 to infinity. I would expect the random number generator to never get close to infinity. Could that cause the consistent bias that you are seeing?
     
  9. May 4, 2016 #8
    Is there any relationship between the two? I mean, the integral is supposed to find the expected value of a random variable numerically. So, in Monte-Carlo simulations, I generated a sufficient number of random samples and then averaged them to find the same quantity.
     
  10. May 4, 2016 #9

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    I guess that depends on how you generate random numbers that can go from 0 to infinity. Are you sure that there is not an upper limit to the numbers that you are generating?
     
  11. May 4, 2016 #10
    The expected value by definition requires an infinity number of samples. But I generated ##10^6## samples, which I think gives close result to infinity, i.e., when I increase from ##10^6## to ##10^7##, the improvement in results is nothing noticeable.
     
  12. May 4, 2016 #11

    nrqed

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    I think that what FactChecker is talking about is the not the number of samples but the range of values sampled by your random generator, i.e. how does one sample fairly all values from 0 to infinity?
     
  13. May 4, 2016 #12
    How can I know that? I use MATLAB for the random sample generator. Since the random variables are exponential random variables with unity mean I use exprnd(1). I'm not sure how it works, and what range of values it generates. But I used it before for other formulas, and it gave accurate results compared to the numerical ones.
     
  14. May 4, 2016 #13

    Dale

    Staff: Mentor

    This is a good thing to check, but it is a check of convergence, not accuracy. Numerical methods are much more subtle than that and require more care if you want high precision.

    It is fundamentally impossible to represent a continuous range of real numbers numerically, due to finite precision. It is further impossible to represent an infinite range.

    Both of these facts will unavoidably lead to a loss of accuracy and precision in any numerical method. So you need to quantify those and determine how to control them. This is well done in NIntegrate, but not in the Monte Carlo simulation. It appears that the Monte Carlo method converges to an inaccurate number.
     
  15. May 5, 2016 #14

    FactChecker

    User Avatar
    Science Advisor
    Gold Member

    That's a pretty good answer. I would expect the MATLAB implementation to be good and you have gotten good results from prior use of it.
     
  16. May 5, 2016 #15
    Then the question still remains: if NIntegrate is accurate, and the implementation of the random number generator is good, what else could be the problem?!
     
  17. May 5, 2016 #16

    ChrisVer

    User Avatar
    Gold Member

    how far away [in standard deviations] are your values?
     
  18. May 5, 2016 #17
    What do you mean? (Excuse me if my question sounded naive)
     
  19. May 5, 2016 #18

    Dale

    Staff: Mentor

    The fact that you have gotten good results from it in other, less demanding, applications does not imply that it is well suited to this application. In particular the numerical errors are not clearly controlled.



    First, you should check the documentation for any known weaknesses. Then generate a large sample and test how much it deviates from an exponential distribution. Then generate a large number of smaller samples to determine if the sampling distribution of the mean is unbiased.
     
  20. May 5, 2016 #19
    Thanks for clarifying. I will test it tomorrow when I get to my office. But if the samples generated represent the exponential distribution fairly well, can we rule out the possibility that the error is due to the random samples generator?
     
  21. May 6, 2016 #20

    ChrisVer

    User Avatar
    Gold Member

    I meant when you calculate an integral with numerical or stochastic [MC] methods, the value is not the exact solution [which may not be known], as a result it comes with an error [itex]\int_a^b f(x) dx = I \pm \delta I[/itex]... I don't know how NIntegrate works and stuff, but for MC you get a statistical error for sure to your integral estimate...
    Also afterall, MC can still give results off since it's a random method... However what I've seen in some cases is that if you take several results out of the MC [itex]\mu_i[/itex] then their average is pretty close to the expected value [within 1 standard deviation]
     
    Last edited: May 6, 2016
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted