2.4 HW 4

1.

This HW in one PDF (letter size)

2.

This HW in one PDF (legal size)

2.4.1 Problem 1

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.

Figure 2.19:problem 1 description

2.4.1.1 Part (a)

\[ f\left ( x\right ) =A\frac{1}{1+x}\] To normalize it, we first 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 find 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 ) $}\]

2.4.1.2 Part(b)

Matlab rand was used to generate number of samples. To see the effect, 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

Algorithm 1:Algorithm to generate random variables from pdf using rand

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

Figure 2.20:Showing how sampling improves with increased \(N\)

Figure 2.21:The case for \(N=10^{6}\) only.

2.4.2 Problem 2

Figure 2.22:problem 2 description

The rejection method was used to first 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 effect. 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

Algorithm 2:Algorithm to generate random variables from pdf using rejection method

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.4.2.1 High LET radiation

Figure 2.23:High LET radiation, \(10^{4}\) trials, normal plot

Figure 2.24:High LET radiation, \(10^{4}\) trials log plot

Figure 2.25:High LET radiation, \(10^{5}\) trials, normal plot

Figure 2.26:High LET radiation, \(10^{5}\) trials log plot

Figure 2.27:High LET radiation, \(10^{6}\) trials, normal plot

Figure 2.28:High LET radiation, \(10^{6}\) trials log plot

Figure 2.29:High LET radiation, \(10^{7}\) trials, normal plot

Figure 2.30:High LET radiation, \(10^{7}\) trials log plot

2.4.2.2 Low LET radiation

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

Figure 2.31:Low LET radiation, \(10^{4}\) trials, normal plot

Figure 2.32:Low LET radiation, \(10^{4}\) trials log plot

Figure 2.33:Low LET radiation, \(10^{5}\) trials, normal plot

Figure 2.34:Low LET radiation, \(10^{5}\) trials log plot

Figure 2.35:Low LET radiation, \(10^{6}\) trials, normal plot

Figure 2.36:Low LET radiation, \(10^{6}\) trials log plot

Figure 2.37:Low LET radiation, \(10^{7}\) trials, normal plot

Figure 2.38:Low LET radiation, \(10^{7}\) trials log plot

2.4.3 Problem 3

Figure 2.39:problem 3 description

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.

Figure 2.40:\(10^{4}\) trials, normal plot

Figure 2.41:\(10^{5}\) trials, normal plot

Figure 2.42:\(10^{6}\) trials, normal plot

Figure 2.43:\(10^{7}\) trials, normal plot

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.

Figure 2.44:\(\{10^{4},10^{5},10^{6},10^{7},10^{8}\}\) trials all on one plot.