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(2pir2))^1/2
(r1, r2: random numbers)}
Any suggestions?
box_muller<-function(n){
RépondreSupprimeru1<-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)
}