mercredi 23 mai 2018

Python, Box-muller Gaussian random number generator & plot his

I am trying to write code in python using the Box-Muller equations, but I do not know how to start!

This is the example I am trying to solve:

{ A peak observed at 900 keV shows a FWHM of 2 keV. Using a Gaussian sampling method listed below, generate 15,000 counts corresponding to the 900 keV peak and save the sampled energies. Create and plot a histogram with a bin width of 0.2 keV and compare with the Gaussian function with the same peak area. Using a data analysis software, try a Gaussian fit to the Monte Carlo data and see the result is close enough to the peak model. Box-Muller Gaussian sampling method: [ Note, the two sampled variable y1, y2 below are for the unit Gaussian distribution (i.e. mu=0,segma=1).

y1 = (-2 ln r1 cos(2pir2))^1/2
y2 = (-2 ln r1 sin(2pi
r2))^1/2

(r1, r2: random numbers)}

Any suggestions?




1 commentaire:

  1. box_muller<-function(n){
    u1<-runif(n)
    r<--2*log(u1)
    u2<-runif(n)
    th<-2*pi*u2
    z1<-sqrt(r)*cos(th)
    z2<-sqrt(r)*sin(th)
    z<-vector(length=n*2)
    i<-1
    j<-1
    while (i<=n*2){
    z[i]=z1[j]
    i=i+1
    z[i]=z2[j]
    j=j+1
    i=i+1
    }
    ratio<-1-(n/j-1)
    print(paste("le ratio de rejet est de ", ratio, "%"))
    return(z)
    }

    RépondreSupprimer