13  Introduction to Numerical Methods in R

Q13.1

The natural logarithm and exponential functions are inverses of each other, so that mathematically \(\log(\exp x) = \exp(\log x) = x\). Show by example that this property does not hold exactly in computer arithmetic. Does the identity hold with near equality? (See all.equal.)

Q13.2

Suppose that \(X\) and \(Y\) are independent random variables, \(X \sim \textrm{Beta}(a,b)\) and \(Y \sim \textrm{Beta}(r,s)\). Then it can be shown that \[ P(X<Y) = \sum_{k=\textrm{max}(r-b,0)}^{r-1} \frac{ {r+s-1 \choose k}{a+b-1 \choose a+r-1-k} }{a+b+r+s-2 \choose a+r-1} \]

Write a function to compute \(P(X<Y)\) for any \(a\),\(b\),\(r\),\(s\) \(>0\). Compare your result with a Monte Carlo estimate of \(P(X<Y)\) for \((a,b) = (10,20)\) and \((r,s) = (5,5)\).

Q13.3

  1. Write a function to compute the \(k^\textrm{th}\) term in \[ \sum_{k=0}^\infty \frac{(-1)^k}{k!2^k} \frac{||a||^{2k+2}}{(2k+1)(2k+2)} \frac{\Gamma(\frac{d+1}{2})\Gamma(k+\frac{3}{2})}{\Gamma(k+\frac{d}{2}+1)} \] where \(d \ge 1\) is an integer, \(a\) is a vector in \(\mathbb{R}^d\), and \(|| \cdot ||\) denotes the Euclidean norm. Perform the arithmetic so that the coefficients can be computed for (almost) arbitrarily large \(k\) and \(d\). (This sum converges for all \(a \in \mathbb{R}^d\).)

  2. Modify the function so that it computes and returns the sum.

  3. Evaluate the sum when \(a=(1,2)^T\).

Q13.4

Find the intersection points \(A(k)\) in \((0,\sqrt{K})\) of the curves \[ S_{k-1}(a) = P \left( t(k-1) > \sqrt{\frac{a^2(k-1)}{k-a^2}} \right) \] and \[ S_{k}(a) = P \left( t(k) > \sqrt{\frac{a^2k}{k+1-a^2}} \right) \] for \(k=4:25\), \(100\), \(500\), \(1000\), where \(t(k)\) is a Student \(t\) random variable with \(k\) degrees of freedom. (These intersection points determine the critical values for a \(t\)-test for scale-mixture errors proposed by Szekely.)

Q13.5

Write a function to solve the equation \[ \frac{2\Gamma(\frac{k}{2})}{\sqrt{\pi(k-1)}\Gamma(\frac{k-1}{2})} \int_0^{c_{k-1}} \left( 1 + \frac{u^2}{k-1} \right)^{-k/2} du = \frac{2\Gamma(\frac{k+1}{2})}{\sqrt{\pi k}\Gamma(\frac{k}{2})} \int_0^{c_k} \left( 1 + \frac{u^2}{k} \right)^{(k+1)/2} du \] for \(a\), where \[ c_k = \sqrt{\frac{a^2k}{k+1-a^2}} \]

Compare the solutions with the points \(A(k)\) in Exercise 13.4.

Q13.6

Write a function to compute the cdf of the Cauchy distribution, which has density \[ \frac{1}{\theta\pi(1+[(x-\eta)/\theta]^2)}, \hspace{2em} -\infty < x < \infty \] where \(\theta>0\). Compare your results to the results from the R function pcauchy. (Also see the source code in pcauchy.c.)