Probability distribution confusion

[yuiop]

Hi, I am trying to write a Monte Carlo type simulation. I have a function for the probability (p) of a photon having an angle (a) between \theta and \pi/4 defined by:

p = \int_{a= \theta}^{a=\pi/2} \cos(2a),\qquad da

...where \theta will normally have a value between -\pi/4 and \pi/4. One immediate problem with the above equation is that the probability of producing a photon of a given exact angle is exactly zero. Do I need to differentiate the equation? Now what I want to do is use a random number generator to produce photons with angles that will match the distribution above when large numbers of photons are produced, but I am having difficulties achieving this. Can anyone offer any suggestions to help me out here?

P.S. If I have posted this question in the wrong forum please advise. I was wondering if it should be in the calculus forum?

 

[EnumaElish]

You need two things:

A. a uniform pseudo-random generator (I'll call this uniform random variable u)

B. a closed-form solution for your integral, F[θ], that you can invert at least locally (i.e. for -π/4 < θ < π/4).

I will assume that you have both A and B, and that F is a well-behaved CDF (i.e. monotone increasing, F[-π/4] = 0, F[π/4] = 1) -- you should verify these. Then:

1. generate a list of draws from u, {u1, ..., un}

2. for each u value, calculate θj = F^-1[uj] to obtain a list of θ values. Those values will be distributed F.

If you don't have a closed-form integral then you may have to numerically calculate each θj.

 

[yuiop]

You need two things:

A. a uniform pseudo-random generator (I'll call this uniform random variable u)

Got one of those

B. a closed-form solution for your integral, F[θ], that you can invert at least locally (i.e. for -π/4 < θ < π/4).

Not quite sure what you mean here. The indefinite integral of cos(2*a) is sin(a)*cos(a). Not quite sure where to go from there either. If I call it u = sin(a)*cos(a) and try to solve for a, I get a fairly complicated result http://www.wolframalpha.com/input/?i=solve+u=+sin(a)+cos(a)+1/2+for+a. Is that what you are referring to?

I will assume that you have both A and B, and that F is a well-behaved CDF (i.e. monotone increasing, F[-π/4] = 0, F[π/4] = 1) -- you should verify these.

I have plot for u = sin(a)*cos(a)+1/2 that smoothly increases from 0 to 1 over that range. It can be seen here
http://www.wolframalpha.com/input/?i=plot+cos(2*a),+sin(a)+cos(a)+1/2+from+a=-pi/4+to+pi/4 alongside a plot for cos(2*a) which does not match your criteria.

Then:

1. generate a list of draws from u, {u1, ..., un}

2. for each u value, calculate θj = F^-1[uj] to obtain a list of θ values. Those values will be distributed F.

If you don't have a closed-form integral then you may have to numerically calculate each θj.

Is the complicated solution for a in the Wolfram Alpha link, the F^-1[uj] you speak of?

 

[yuiop]

You need two things:

A. a uniform pseudo-random generator (I'll call this uniform random variable u)

B. a closed-form solution for your integral, F[θ], that you can invert at least locally (i.e. for -π/4 < θ < π/4).

I will assume that you have both A and B, and that F is a well-behaved CDF (i.e. monotone increasing, F[-π/4] = 0, F[π/4] = 1) -- you should verify these. Then:

1. generate a list of draws from u, {u1, ..., un}

2. for each u value, calculate θj = F^-1[uj] to obtain a list of θ values. Those values will be distributed F.

If you don't have a closed-form integral then you may have to numerically calculate each θj.

O.K. I managed to get a slightly simpler solution http://www.wolframalpha.com/input/?i=solve+v=+sin(a)+cos(a)+for+a by using v=cos(a)*sin(a) and solving for a, using v as a substitution for u-1/2. This produced a good match to the theoretical distribution using pseudo-random data. Thanks for the help. You obviously know what you are doing. Cheers!