- 1.
- This HW in one PDF (letter size)
- 2.
- This HW in one PDF (legal size)

Important note on Matlab: I used Matlab 2016a to write all the software. In particular, the function histcounts() was used for binning as recommended by Matlab help pages. This function do not exist in Matlab 2014a and was added in Matlab 2014b.

\[ f\left ( x\right ) =A\frac{1}{1+x}\] To normalize it, we ﬁrst solve for \(A\) from\[ \int _{-\infty }^{\infty }f\left ( x\right ) dx=1 \] Hence\begin{align*} 1 & =\int _{-\infty }^{\infty }A\frac{1}{1+x}dx\\ & =A\left [ \arctan \left ( x\right ) \right ] _{-\infty }^{\infty }\\ & =A\left ( \frac{\pi }{2}+\frac{\pi }{2}\right ) \end{align*}

Therefore \[ \fbox{$A=\frac{1}{\pi }$}\] Now we need to ﬁnd cumulative probability distribution \(F\left ( x\right ) \) and invert then it. \begin{align*} F\left ( x\right ) & =\int _{-\infty }^{x}f\left ( x\right ) dx\\ & =\int _{-\infty }^{x}\frac{1}{\pi }\frac{1}{1+x}dx\\ & =\frac{1}{\pi }\left ( \arctan \left ( x\right ) \right ) _{-\infty }^{x}\\ & =\frac{1}{\pi }\left ( \arctan \left ( x\right ) +\frac{\pi }{2}\right ) \end{align*}

Hence\[ F\left ( x\right ) =\frac{1}{\pi }\arctan \left ( x\right ) +\frac{1}{2}\] To invert, \[ \pi \left ( F\left ( x\right ) -\frac{1}{2}\right ) =\arctan \left ( x\right ) \] Hence \[ \fbox{$x=F^-1\left ( x\right ) =\tan \left ( \pi \left ( F\left ( x\right ) -\frac{1}{2}\right ) \right ) $}\]

Matlab rand was used to generate number of samples. To see the eﬀect, the sampling was increased from \(10^{2}\) to \(10^{6}\). The area under the pdf generated using sampling was normalized so that its area is one. The following plots show the result. The source code is included.

The algorithm used is the following

1:
procedure generate_RV

2: \(F(x) \gets \text{cumulative distribution}\)

3: \(N \gets \text{Number of trials}\)

4: for \(i \gets 1,N\) do

5: \(x(i) \gets \text{rand()}\)

6: \(data(i) \gets F^{-1}(x(i))\)

7: end for

8: generate normalized histogram from \(data\)

9: end procedure

2: \(F(x) \gets \text{cumulative distribution}\)

3: \(N \gets \text{Number of trials}\)

4: for \(i \gets 1,N\) do

5: \(x(i) \gets \text{rand()}\)

6: \(data(i) \gets F^{-1}(x(i))\)

7: end for

8: generate normalized histogram from \(data\)

9: end procedure

The rejection method was used to ﬁrst generate \(f_{H}(D)\) and then was used again to generate \(F_{L}(D)\). So the result shows each of these separately below.

For each case, number of trials was increased to see the eﬀect. The following trials were used \(\{10^{4},10^{5},10^{6},10^{7}\}\). The last amount is the one asked for in this problem. For each case, both the similogy plot is shown and the normal scale plot is also shown. These results shows the rejection method improves as more trials are used. The analytic pdf is also plotted against the generated one to compare.

Matlab source code is included. Implemented on Matlab 2016a.

The algorithm used is the following

1:
procedure generate_RV

2: \(f(x) \gets \text{probability distribution}\)

3: \(N \gets \text{Number of trials}\)

4: \(k \gets 0\)

5: for \(i \gets 1,N\) do

6: \(\zeta _1 \gets \text{rand()}\)

7: \(\zeta _2 \gets \text{rand()}\)

8: if \(\zeta _2 \leq f(\zeta _1)\) then

9: \(k\gets k+1\)

10: \(\text{count[k]} \gets \zeta _1\)

11: end if

12: end for

13: generate normalized histogram from count

14: end procedure

2: \(f(x) \gets \text{probability distribution}\)

3: \(N \gets \text{Number of trials}\)

4: \(k \gets 0\)

5: for \(i \gets 1,N\) do

6: \(\zeta _1 \gets \text{rand()}\)

7: \(\zeta _2 \gets \text{rand()}\)

8: if \(\zeta _2 \leq f(\zeta _1)\) then

9: \(k\gets k+1\)

10: \(\text{count[k]} \gets \zeta _1\)

11: end if

12: end for

13: generate normalized histogram from count

14: end procedure

The above was repeated for the low LET radiation. The result is given below

The following trials were used \(\{10^{4},10^{5},10^{6},10^{7}\}\). The last amount is the one asked for in this problem. For each case, analytic pdf is plotted against the generated one to compare. As more trials used, the approximation to the true pdf is improved.

The following plots adds \(10^{8}\) number of trials in order to see how close the generated pdf will be to the theoretical pdf in the plateau region (the region between 15 and 45 years). Generated pdf becomes closer to the theoretical pdf, but there always be a small gap at the top. With more trials, the generated pdf will converge to the theoretical pdf.