Mathematica Numerical vs. Monte-Carlo Simulations

[EngWiPy]

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?

 

[EngWiPy]

How can I upload an image from my computer?

 

[Dale]

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.

 

[EngWiPy]

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.

 

[EngWiPy]

See here, please.

 

[Dale]

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.

 

[FactChecker]

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?

 

[EngWiPy]

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?

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.

 

[FactChecker]

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.

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?

 

[EngWiPy]

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?

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.

 

[nrqed]

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.

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?

 

[EngWiPy]

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?

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.

 

[Dale]

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.

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.

 

[FactChecker]

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.

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.

 

[EngWiPy]

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?!

 

[ChrisVer]

how far away [in standard deviations] are your values?

 

[EngWiPy]

how far away [in standard deviations] are your values?

What do you mean? (Excuse me if my question sounded naive)

 

[Dale]

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?!

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.

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.

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.

 

[EngWiPy]

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?

 

[ChrisVer]

What do you mean? (Excuse me if my question sounded naive)

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 \int_a^b f(x) dx = I \pm \delta I... 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 \mu_i then their average is pretty close to the expected value [within 1 standard deviation]

 

[OrangeDog]

I think your monte carlo equation is screwed up. You can't get the solution you are seeking unless the exponential equals 0, which can't happen without an imaginary number.

 

[Hepth]

Test NIntegrate with different methods

NIntegrate[___, Method-> "MonteCarlo"]
"GlobalAdaptive" [default]
"DuffyCoordinates"
"MonteCarlo"
"QuasiMonteCarlo"
"AdaptiveMonteCarlo"
"AdaptiveQuasiMonteCarlo"
"DoubleExponential"

you might have to add MaxPoints-> 10^7

Plain montecarlo seems to do the worst.

 

[EngWiPy]

I think your monte carlo equation is screwed up. You can't get the solution you are seeking unless the exponential equals 0, which can't happen without an imaginary number.

I'm sorry. I didn't understand what you mean. Again, why my Monte Carlo simulations is screwed up?

 

[EngWiPy]

Test NIntegrate with different methods

NIntegrate[___, Method-> "MonteCarlo"]
"GlobalAdaptive" [default]
"DuffyCoordinates"
"MonteCarlo"
"QuasiMonteCarlo"
"AdaptiveMonteCarlo"
"AdaptiveQuasiMonteCarlo"
"DoubleExponential"

you might have to add MaxPoints-> 10^7

Plain montecarlo seems to do the worst.

I'm doing Monte Carlo simulations on MATLAB. What does the above code do?

 

[EngWiPy]

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.

I tried to check the generator. But after I draw the exponential distribution for an exponential random variable with mean 1, I didn't know how to generate the exponential random variables and match them to the numerical result. I wrote the following:

x=0:0.1:5;
y=exp(-x);

ySim=exprnd(1,1,length(x));%This generates an 1-by-length(x) array of exponential random variables of mean 1

plot(x,y,x,ySim);

I got the results attached. I think I didn't do it right, did I?

 

[ChrisVer]

why don't you create some random numbers in the range (a,b) and then pass them in the function f(x) you want ( like the exponential )?

 

[EngWiPy]

why don't you create some random numbers in the range (a,b) and then pass them in the function f(x) you want ( like the exponential )?

What would the range be in my case? and how will this verify the exprnd generator which generates exponential random variables I use in my simulations? I don't generate random variables and pass them to an exponential function.

 

[Dale]

I tried to check the generator. But after I draw the exponential distribution for an exponential random variable with mean 1, I didn't know how to generate the exponential random variables and match them to the numerical result. I wrote the following:

x=0:0.1:5;
y=exp(-x);

ySim=exprnd(1,1,length(x));%This generates an 1-by-length(x) array of exponential random variables of mean 1

plot(x,y,x,ySim);

I got the results attached. I think I didn't do it right, did I?

You will want to plot the random data using either a histogram or a qq plot.

The histogram will probably look pretty good, but the qq plot will be more sensitive to deviations.

 

[EngWiPy]

You will want to plot the random data using either a histogram or a qq plot.

The histogram will probably look pretty good, but the qq plot will be more sensitive to deviations.

How can I do it with the histogram? I have the exponential random variables in ySim.

 

[Dale]

I don't know Matlab, I use Mathematica and R. I would think that Matlab has commands for plotting histograms and qq plots, but I cannot tell you their syntax.

 

[EngWiPy]

It's histogram(x). When I used with the data I have, I got the attached figure. Does that make sense?

 

[Dale]

It's histogram(x). When I used with the data I have, I got the attached figure. Does that make sense?

Yes, that seems right. Notice how it does not go out to infinity and how it is not strictly monotonically decreasing. The variations are small, but they will make a difference in achieving a given level of numerical precision.

 

[EngWiPy]

Yes, that seems right. Notice how it does not go out to infinity and how it is not strictly monotonically decreasing. The variations are small, but they will make a difference in achieving a given level of numerical precision.

Can I do anything about it? I'm writing a scientific journal, do you think such justification is acceptable for why there is a difference between the two curves?

 

[Dale]

Hmm, OK, if you are writing a scientific paper then you want to hold your work to a higher standard. Let's look at a couple more things.

First, instead of generating 1,000,000 numbers that are exponentially distributed and plotting the histogram of the distribution, let's investigate the sampling distribution. Instead of 1,000,000 numbers, generate 1,000 sets of 1,000 numbers each. For each set calculate the mean and the standard deviation (1,000 means and 1,000 standard deviations) and then plot the histograms of those. Let's see if there is any bias.

Second, see if you can easily get a qq plot in Matlab.

 

[EngWiPy]

In the first figure, I didn't generate $10^6$ samples. In the attached figure I did. Does this change anything before going to implement your suggestion?

 

[Dale]

In the first figure, I didn't generate $10^6$ samples. In the attached figure I did. Does this change anything before going to implement your suggestion?

No. It still shows the same key features, but it is smaller and harder to see. The range is not infinite and the curve is not smooth. This is to be expected in any random number generator, but it is a source of numerical error.

If we are lucky then the sampling distribution will show some small bias.

 

[EngWiPy]

OK, I did your suggestion, and got the attached figures for the mean and standard deviation. What can we infer from this?

 

[EngWiPy]

Also attached are the qq plot for the means and standard deviations calculated above

 

[Dale]

So it looks like the RNG is not too bad but the sampling distribution does show a bit heavier tails than normal and a bit of right skew.

Can you do a qq plot of the ySim against the exponential distribution instead of the normal distribution?

Also, what is the mean and standard deviation of the sampling distribution of the mean?

 

[EngWiPy]

So it looks like the RNG is not too bad but the sampling distribution does show a bit heavier tails than normal and a bit of right skew.

Q1: Can you do a qq plot of the ySim against the exponential distribution instead of the normal distribution?

Q2: Also, what is the mean and standard deviation of the sampling distribution of the mean?

I didn't understand your questions. For Q1: What do you mean by qq plotting ySim against the exponential distribution? and how? I qq plot ySim alone and got the attached figure.

I also didn't understand Q2. Sorry, but can you elaborate, please.

 

[EngWiPy]

OK, for Q2: the mean and the standard deviation of the mean are 1.0001 and 0.0309, respectively.

 

[Dale]

OK, for Q2: the mean and the standard deviation of the mean are 1.0001 and 0.0309, respectively.

OK, that is essentially unbiased.

 

[EngWiPy]

So, what is the implication of this on my result?

 

[Dale]

I didn't understand your questions. For Q1: What do you mean by qq plotting ySim against the exponential distribution? and how? I qq plot ySim alone and got the attached figure

So the qq plots you have posted so far plot against the normal distribution. The ySim values are not normally distributed so the qq plot should be done against the exponential distribution.

 

[Dale]

It should be some sort of an option to plot against different distributions.

 

[EngWiPy]

Right. I did it, and got the attached figure.

 

[Dale]

Good. Do you see the deviations on the upper right and how there is some discrete steps where there are 4 points with the same vertical position? The errors from that are not well characterized. So I think you can justify not using the Monte Carlo.

It appears that there is not much, if any, bias in the sampling distribution, but the asymmetric numerical errors make those results suspicious.

 

[EngWiPy]

I need to use Monte Carlo simulations as an independent verification of my results. To the best of my knowledge in my field, unless there is some numerical approximation in the mathematical analysis, Monte Carlo simulations give accurate results. I'm not sure why it's not the case here. I'm concerned that my paper will be rejected because of this simulation only, as I think the mathematical analysis is correct. I forgot if we discussed this before, but is the Monte Carlo simulation procedure I described in the other thread (see post #5) correct?

 

[EngWiPy]

Any comment on this? Is it acceptable to say on a scientific article that the small discrepancy in results is due to some inaccuracy in the simulations?

 

[mfig]

Is it acceptable to say on a scientific article that the small discrepancy in results is due to some inaccuracy in the simulations?

What is the small discrepancy? 1 part in 1 million? 1 part in 1 billion? Is there some reason you are using an exponential distribution and not a uniform distribution for your sampling? If you post just the MC integration code, we might be able to see a problem with it.