Monte carlo algorithms in C++

  • #1
Chromatic_Universe
Gold Member
13
0
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.
 

Answers and Replies

  • #2
Nguyen Son
54
36
Could you code Matlab, besides C++?
 
  • #3
Chromatic_Universe
Gold Member
13
0
Could you code Matlab, besides C++?
Used it once before, but lost touch.
 
  • #4
Hello @Chromatic_Universe
:welcome:
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.
 
  • #5
Chromatic_Universe
Gold Member
13
0
Hello @Chromatic_Universe
:welcome:

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 !
 
  • #6
Chromatic_Universe
Gold Member
13
0
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.
 
  • #7
FactChecker
Science Advisor
Homework Helper
Gold Member
7,593
3,314
We have the function sin^2/x^2, given x>=0.
Is this correct? What are you taking the sin of? sin2(?)/x2
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.)
 
  • #8
Chromatic_Universe
Gold Member
13
0
Is this correct? What are you taking the sin of? sin2(?)/x2
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)2/x2. 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)2/x2(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.
 
  • #9
FactChecker
Science Advisor
Homework Helper
Gold Member
7,593
3,314
@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)2/x2(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.
 
  • #10
FactChecker
Science Advisor
Homework Helper
Gold Member
7,593
3,314
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 x0, where P( x>x0) = 1/100. Then use the rejection method for 0<x<x0 and randomly (with probability 1/100) use the rejection method on 100*f(x) where x > x0.
You can extend this to more strata if necessary.
 
Last edited:

Suggested for: Monte carlo algorithms in C++

Replies
1
Views
2K
  • Last Post
Replies
2
Views
622
Replies
9
Views
955
Replies
11
Views
2K
Replies
20
Views
1K
Replies
1
Views
517
Replies
18
Views
731
Replies
2
Views
491
Top