Need Help with Monte Carlo Algorithms in C++?

[Chromatic_Universe]

I would like to discuss code for hit and miss monte carlo methods, and also monte carlo with veto algorithm in C++. Since I am coding in C++ after a long time, I am messed up with syntax too. I have a specific set of problems to work with. If interested we can start working on it here.

 

[Nguyen Son]

Could you code Matlab, besides C++?

 

[Chromatic_Universe]

Could you code Matlab, besides C++?

Used it once before, but lost touch.

 

[NFuller]

Hello @Chromatic_Universe

I have a specific set of problems to work with.

There are numerous members here (including myself) with experience in C++ and monte carlo methods. If you have specific problems that you are stuck on, please create a post for each specific set of issues. Also if these are homework problems, please remember to use the homework template in your post.

 

[Chromatic_Universe]

Hello @Chromatic_Universe


There are numerous members here (including myself) with experience in C++ and monte carlo methods. If you have specific problems that you are stuck on, please create a post for each specific set of issues. Also if these are homework problems, please remember to use the homework template in your post.

Will do so! Thanks !

 

[Chromatic_Universe]

We have the function sin^2/x^2, given x>=0. Now we have to generate random numbers according to the given distribution using Monte Carlo hit-and-miss method. I want to do it in C++.
@NFuller you can see the problem now.

 

[FactChecker]

We have the function sin^2/x^2, given x>=0.

Is this correct? What are you taking the sin of? sin²(?)/x²

Now we have to generate random numbers according to the given distribution using Monte Carlo hit-and-miss method. I want to do it in C++.

Do you mean the Monte Carlo rejection sampling method for generating a random variable from a probability density? This example presents a problem because of the unbounded x range. That is more than a programming problem. Do you have a specific way that you want to address that problem? (You may want to consider that there will be a lot of samples rejected when x is large.)

 

[Chromatic_Universe]

Is this correct? What are you taking the sin of? sin²(?)/x²
Do you mean the Monte Carlo rejection sampling method for generating a random variable from a probability density? This example presents a problem because of the unbounded x range. That is more than a programming problem. Do you have a specific way that you want to address that problem? (You may want to consider that there will be a lot of samples rejected when x is large.)

@FactChecker It would be (sin x)²/x². And it does not blow up at x=0 because we apply L'Hopital's rule(https://en.wikipedia.org/wiki/L'Hôpital's_rule), hence it's finite(You can also check the intensity distribution for a single slit diffraction experiment - http://hyperphysics.phy-astr.gsu.edu/hbase/phyopt/sinint.html).
A way to do it would be to take a test function g(y) which is always greater than f(x)=(sin x)²/x²(where x is a random number) for the given range and generate random "y' " and calculate values of g(y') and check with f(x). If g(y) <= f(x), we accept, otherwise we reject it. Then we repeat the procedure again.
We consider g(y) to increase efficiency of the algorithm. Otherwise, as you mentioned -

You may want to consider that there will be a lot of samples rejected when x is large.

 

[FactChecker]

@FactChecker It would be (sin x)^2/x^2. And it does not blow up at x=0 because we apply L'Hopital's rule(https://en.wikipedia.org/wiki/L'Hôpital's_rule), hence it's finite(You can also check the intensity distribution for a single slit diffraction experiment - http://hyperphysics.phy-astr.gsu.edu/hbase/phyopt/sinint.html).

Yes. I forgot for a moment that the sin was squared. I corrected it a few minutes later.

A way to do it would be to take a test function g(y) which is always greater than f(x)=(sin x)²/x²(where x is a random number) for the given range and generate random "y' " and calculate values of g(y') and check with f(x). If g(y) <= f(x), we accept, otherwise we reject it. Then we repeat the procedure again.
We consider g(y) to increase efficiency of the algorithm. Otherwise, as you mentioned -

Ok. You may want to consider the problems presented by the unbounded x range and the high number of rejections for large values of x.

 

[FactChecker]

I don't see at the moment how you can use g(x) without changing the probability density.

You can try a simple rejection method using reasonably high limit on x and see if there is a problem with too many rejections.
If the rejections are too high, I suggest a "stratified" rejection method. For instance, figure out x₀, where P( x>x₀) = 1/100. Then use the rejection method for 0<x<x₀ and randomly (with probability 1/100) use the rejection method on 100*f(x) where x > x₀.
You can extend this to more strata if necessary.

 

[NFuller]

It may be helpful to take a look here.
http://home.fnal.gov/~mrenna/lu_update/lutp0613man2/node40.html
Is this what you are trying to do?