Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
PDF to CDF in MATLAB

Subject: PDF to CDF in MATLAB

From: Hemming

Date: 13 Jan, 2013 19:40:09

Message: 1 of 7

Hi!
Im trying to extract a scattering angle for a photon using the Klein-Nishina scattering angle distribution (KN in the code) and for this I need the CDF (of KN) to be able to use the Monte Carlo method when that is achieved. All i've managed so far is to plot the PDF between 0 degrees and Pi to see that it looks alright, and that it has that "peanut shape". Ive tried to use the built in CDF function but it seems very slow.

E_gamma=0.140;
alpha=0.511;
P=zeros(3142,1);
KN_matrix=zeros(3142, 1);

for k=1:3142
j=(k-1)/1000;
PE=1/(1+((E_gamma/alpha)*(1-cos(j))));
KN=(PE^2*(PE+(1/PE)-1+(cos(j))^2))/2;
KN_matrix(k,1)=(KN);
plot(j,KN_matrix(k,1))
hold on
end

Thanks!
axis equal

Subject: PDF to CDF in MATLAB

From: TideMan

Date: 13 Jan, 2013 20:55:11

Message: 2 of 7

On Monday, January 14, 2013 8:40:09 AM UTC+13, Hemming wrote:
> Hi!
>
> Im trying to extract a scattering angle for a photon using the Klein-Nishina scattering angle distribution (KN in the code) and for this I need the CDF (of KN) to be able to use the Monte Carlo method when that is achieved. All i've managed so far is to plot the PDF between 0 degrees and Pi to see that it looks alright, and that it has that "peanut shape". Ive tried to use the built in CDF function but it seems very slow.
>
>
>
> E_gamma=0.140;
>
> alpha=0.511;
>
> P=zeros(3142,1);
>
> KN_matrix=zeros(3142, 1);
>
>
>
> for k=1:3142
>
> j=(k-1)/1000;
>
> PE=1/(1+((E_gamma/alpha)*(1-cos(j))));
>
> KN=(PE^2*(PE+(1/PE)-1+(cos(j))^2))/2;
>
> KN_matrix(k,1)=(KN);
>
> plot(j,KN_matrix(k,1))
>
> hold on
>
> end
>
>
>
> Thanks!
>
> axis equal

Instead of plotting in a loop, define vectors:
j=[0:3141]'/1000;
PE=1/(1+((E_gamma/alpha)*(1-cos(j))));
KN_matrix=(PE.^2.*(PE+(1./PE)-1+(cos(j)).^2))/2;
plot(j,KN_matrix)

Notice the dots before the operators to get element-by-element multiplication, exponentiation, etc.

Subject: PDF to CDF in MATLAB

From: Roger Stafford

Date: 13 Jan, 2013 22:02:11

Message: 3 of 7

"Hemming" wrote in message <kcv2ip$mv1$1@newscl01ah.mathworks.com>...
> Im trying to extract a scattering angle for a photon using the Klein-Nishina scattering angle distribution (KN in the code) and for this I need the CDF (of KN) to be able to use the Monte Carlo method when that is achieved. All i've managed so far is to plot the PDF between 0 degrees and Pi to see that it looks alright, and that it has that "peanut shape". Ive tried to use the built in CDF function but it seems very slow.
- - - - - - - - - -
  Do I understand correctly that this is a probability density function with respect to a solid angle of scattering? Wouldn't that need to be known in obtaining a cumulative distribution? Perhaps you need to explain your problem in much greater detail since you are communicating with people who are not necessarily conversant with quantum electrodynamics. Pretend we know nothing about the physics of the subject and explain your question in purely mathematical terms and in great detail.

Roger Stafford

Subject: PDF to CDF in MATLAB

From: Roger Stafford

Date: 14 Jan, 2013 01:24:12

Message: 4 of 7

"Roger Stafford" wrote in message <kcvat3$k7t$1@newscl01ah.mathworks.com>...
> "Hemming" wrote in message <kcv2ip$mv1$1@newscl01ah.mathworks.com>...
> > Im trying to extract a scattering angle for a photon using the Klein-Nishina scattering angle distribution (KN in the code) and for this I need the CDF (of KN) to be able to use the Monte Carlo method when that is achieved. All i've managed so far is to plot the PDF between 0 degrees and Pi to see that it looks alright, and that it has that "peanut shape". Ive tried to use the built in CDF function but it seems very slow.
> - - - - - - - - - -
> Do I understand correctly that this is a probability density function with respect to a solid angle of scattering? Wouldn't that need to be known in obtaining a cumulative distribution? .......
- - - - - - - - - -
  One more thought about your question. If it is solid angle that the probability density is taken with respect to, then its cumulative distribution in terms of scatter angle is easy to find. However if it is this scattering you wish to simulate in a Monte Carlo process using matlab's 'rand' generator, you would need the inverse of this CDF function and that may not be so easy to calculate. There are some functions in matlab that may able to help in this. There is also the alternative of a random procedure using rejection to achieve the proper simulation distribution. Anyway please give us the details of what it is you wish to do.

Roger Stafford

Subject: PDF to CDF in MATLAB

From: Hemming

Date: 14 Jan, 2013 07:54:08

Message: 5 of 7

Thanks for all the replies, i realize my formulation was very vague. Here's another go:

My problem is that different scattering angles of Compton scattering have different propabilities. High energy photons is, for instance, more likely to scatter 2 degrees than 90. This is were Klein-Nishina (KN) comes into play. I would like to create a distribution (a CDF), using the KN-function between 0-180 degrees (which is my PDF).
So what i should do is integrate KN with respect to the scattering angle (j) (create the inverse of this function as you put it Roger) and this has proved difficult by hand.
If this could be achieved i would use this CDF to randomly select points in the distribution representing different angles.

I hope this clears it up a bit, otherwise i'll give it another go. Thanks again!

Hemming


"Roger Stafford" wrote in message <kcvmns$rgv$1@newscl01ah.mathworks.com>...
> "Roger Stafford" wrote in message <kcvat3$k7t$1@newscl01ah.mathworks.com>...
> > "Hemming" wrote in message <kcv2ip$mv1$1@newscl01ah.mathworks.com>...
> > > Im trying to extract a scattering angle for a photon using the Klein-Nishina scattering angle distribution (KN in the code) and for this I need the CDF (of KN) to be able to use the Monte Carlo method when that is achieved. All i've managed so far is to plot the PDF between 0 degrees and Pi to see that it looks alright, and that it has that "peanut shape". Ive tried to use the built in CDF function but it seems very slow.
> > - - - - - - - - - -
> > Do I understand correctly that this is a probability density function with respect to a solid angle of scattering? Wouldn't that need to be known in obtaining a cumulative distribution? .......
> - - - - - - - - - -
> One more thought about your question. If it is solid angle that the probability density is taken with respect to, then its cumulative distribution in terms of scatter angle is easy to find. However if it is this scattering you wish to simulate in a Monte Carlo process using matlab's 'rand' generator, you would need the inverse of this CDF function and that may not be so easy to calculate. There are some functions in matlab that may able to help in this. There is also the alternative of a random procedure using rejection to achieve the proper simulation distribution. Anyway please give us the details of what it is you wish to do.
>
> Roger Stafford

Subject: PDF to CDF in MATLAB

From: Hemming

Date: 14 Jan, 2013 09:52:06

Message: 6 of 7

Hello again!

I just found a 400 page book on Kahn's method in our university library that should propably help me here. I think i should give it a read before troubling you all with my non existent MATLAB skills again. Thanks for time and effort!


"Hemming" wrote in message <kd0dj0$bsu$1@newscl01ah.mathworks.com>...
> Thanks for all the replies, i realize my formulation was very vague. Here's another go:
>
> My problem is that different scattering angles of Compton scattering have different propabilities. High energy photons is, for instance, more likely to scatter 2 degrees than 90. This is were Klein-Nishina (KN) comes into play. I would like to create a distribution (a CDF), using the KN-function between 0-180 degrees (which is my PDF).
> So what i should do is integrate KN with respect to the scattering angle (j) (create the inverse of this function as you put it Roger) and this has proved difficult by hand.
> If this could be achieved i would use this CDF to randomly select points in the distribution representing different angles.
>
> I hope this clears it up a bit, otherwise i'll give it another go. Thanks again!
>
> Hemming
>
>
> "Roger Stafford" wrote in message <kcvmns$rgv$1@newscl01ah.mathworks.com>...
> > "Roger Stafford" wrote in message <kcvat3$k7t$1@newscl01ah.mathworks.com>...
> > > "Hemming" wrote in message <kcv2ip$mv1$1@newscl01ah.mathworks.com>...
> > > > Im trying to extract a scattering angle for a photon using the Klein-Nishina scattering angle distribution (KN in the code) and for this I need the CDF (of KN) to be able to use the Monte Carlo method when that is achieved. All i've managed so far is to plot the PDF between 0 degrees and Pi to see that it looks alright, and that it has that "peanut shape". Ive tried to use the built in CDF function but it seems very slow.
> > > - - - - - - - - - -
> > > Do I understand correctly that this is a probability density function with respect to a solid angle of scattering? Wouldn't that need to be known in obtaining a cumulative distribution? .......
> > - - - - - - - - - -
> > One more thought about your question. If it is solid angle that the probability density is taken with respect to, then its cumulative distribution in terms of scatter angle is easy to find. However if it is this scattering you wish to simulate in a Monte Carlo process using matlab's 'rand' generator, you would need the inverse of this CDF function and that may not be so easy to calculate. There are some functions in matlab that may able to help in this. There is also the alternative of a random procedure using rejection to achieve the proper simulation distribution. Anyway please give us the details of what it is you wish to do.
> >
> > Roger Stafford

Subject: PDF to CDF in MATLAB

From: Roger Stafford

Date: 14 Jan, 2013 09:59:06

Message: 7 of 7

"Hemming" wrote in message <kd0dj0$bsu$1@newscl01ah.mathworks.com>...
> Thanks for all the replies, i realize my formulation was very vague. Here's another go:
>
> My problem is that different scattering angles of Compton scattering have different propabilities. High energy photons is, for instance, more likely to scatter 2 degrees than 90. This is were Klein-Nishina (KN) comes into play. I would like to create a distribution (a CDF), using the KN-function between 0-180 degrees (which is my PDF).
> So what i should do is integrate KN with respect to the scattering angle (j) (create the inverse of this function as you put it Roger) and this has proved difficult by hand.
> If this could be achieved i would use this CDF to randomly select points in the distribution representing different angles.
>
> I hope this clears it up a bit, otherwise i'll give it another go. Thanks again!
>
> Hemming
- - - - - - - - - -
  My understanding is that the Klein-Nishina formula gives a differential cross section for photons which amounts to a relative probability density with respect to solid angle. That is, roughly speaking it is the relative probability of a scattering within some tiny solid angle divided by the amount of that solid angle. This means that to get a cumulative probability distribution it is necessary to take the integral of this Klein-Nishina expression with respect to cos(theta) where theta is the angle of scattering because of the nature of solid angle measure. This is different from integrating with respect to just theta. It also means that it is fairly easy to obtain this integral since if we substitute x = cos(theta), it would then be the integral of a rational function of x with respect to x. This is just a problem in partial fractions which is solved using elementary calculus.
Some of the resulting terms will give rise to other rational functions of x, but one of them yields a logarithm expression. With the proper scaling and substituting cos(theta) again in place of x you would then have a valid cumulative distribution function as a function of theta. As I say, this is fairly easy to do. It does not require any numerical computation.

  However, to do a Monte Carlo simulation of this process using a random number generator which has a uniform distribution like matlab's 'rand' function, it is necessary to find the inverse of the above CDF. By this I mean that the random number generator generates a value and you must then find a value of theta that would give the CDF function that random value. This is the reverse of directly evaluating the CDF given a value of theta. It is the combination of the rational terms and logarithm terms in the CDF expression that makes it unlikely that an explicit formula can be found to do this. Matlab does have functions like 'fzero' which can solve implicit equations but these are somewhat time-consuming and depend on making good initial estimates.

  There are also what are known as 'rejection' methods in which a certain fraction of random numbers generated must be rejected in such a way as to give the appropriate distribution. I would say that as things stand at the moment it looks as though such a rejection method would be your best bet for doing a Monte Carlo procedure of this scattering, in spite of the need to generate considerably more random numbers.

Roger Stafford

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us