Project 1

The projects should be handed in here as a single pdf file no later than Sunday February 9 at 23:59. You should work in groups of ideally two students. Fill in the name and email of your groups in this google docs spreadsheet or fill in only your own name if you want us to assist you in finding another student to work with.

Use LaTeX (you may want to use knitr in combination with LaTeX, this also works in overleaf) (you can use this template as a starting point) or Rmarkdown (template).

All coding should use R. Use variable names in R that corresponds to the mathematical model. Add comments to your code if needed. Also, format your R code properly, see e.g. this guide. Note that not all R-code chuncks needs to be visible in the report, e.g. code chuncks with complex code used to generate figures should be made invisible (use the echo=FALSE or inclue=FALSE chunk options), instead clearly explain what is shown by each Figure.

The report should be a self-contained concisely written document written in either English or any Scandinavian language. Mathematics should be written using proper mathematical notation and punctuation and all symbols should be explained.

The project reports counts 30% towards the final grade in total.

You are allowed (and encouraged) to discuss the different problems between groups but copying someone elses report or contributing to this is not (see regulations on cheating). Also see regulations on the use of AI-tools. Note in particular that you must enclose an AI-declaration form if you have used AI as an aid when preparing your report. This does by no means imply that using AI is recommended.

On Thursday, February 616 at 16:15-18:00 in R60 (in week 6) each group will present one of the problems below orally in class (or via zoom/teams). This presentations are obligatory but do not count towards the final grade. We will assign the different problems to different groups as soon as you you fill in your groups in the above google spread sheet. Thus, the earlier you do that, the more more time to prepare slides etc. for the oral presentation. Upload your slides here at any time before your presentation.

The presentation should give a brief overview of the solution and emphasize difficult parts of the problem e.g. through graphs and key equations. Each presentation should be about 7 minutes plus 3 minutes for questions.

To prepare slides for the oral presentation, we recommend Beamer.

We will provide guidance on Tuesdays at 10:15-12:00 in G21 in weeks 5 and 6. You can also ask questions in the EdStem forum.

Suppose that \(X_1\) has density \(f(x)=\frac1{\pi\sqrt{x(1-x)}}\) for \(0<x<1\).

a) Implement a function in R simulating \(n\) realization of \(X_1\) from a uniform random variable \(U\sim \operatorname{unif}(0,1)\) using the inversion method and simulate \(n=10000\) realizations. Test that your method works by constructing a suitable sanity check.

Consider the stochastic process \(X_1,X_2,\dots\) defined by the (deterministic) recursion \(X_{t+1}=4X_t(1-X_t)\).

b) For each realization of \(X_1\), compute the associated values of \(X_2,X_3,\dots,X_{10}\). You'll probably want to store the realizations in a \(10000 \times 10\) matrix. How does the distribution of each \(X_t\) appear to compare to that of \(X_1\)? Use Monte-Carlo integration to estimate the variance matrix of \(X_1,X_2,\dots,X_{10}\). Comment on what you see.

c) Try to prove the pattern suggested by what you observe in point b) analytically. Hint: Express each \(X_t\) as functions of \(U\) only and use this to obtain \(\operatorname{Cov}(X_t,X_s)\) via the definition of the covariance. Hint: You will perhaps want to use the double-angle formulae for the cosine function

Suppose that \(X,Y\) are positive random variables with joint cumulative distribution function \[ F_{X,Y}(x,y)=\exp\Big(-\frac1x-\frac1y-\frac1{xy}\Big), \] for \(x>0,y>0\).

a) Find a method for simulating from the joint distribution of \(X\) and \(Y\).

Hints: First find \(F_X(x)=P(X\le x)\) and use this to derive a method for simulating \(X\). Next, let \(F_{Y|X}(y|x)\) denote the conditional cumulative distribution of \(Y\) given \(X=x\). Show that \[ F_{Y|X}(y|x)=\frac1{f_X(x)}\frac{\partial}{\partial x}F_{X,Y}(x,y). \] Use this relation to show that conditional distribution of \(Y\) conditional on \(X=x\) is a mixture of two inverse Gamma distributions (with weights and parameters that depend on \(x\)), that is, \(f_{Y|X}(y|x)\) can be made equal to \(\sum_{k=1}^2 w_k f_k(y;\alpha_k,\beta)\) where \(f_k\) are pdfs (including their normalising constants) of inverse gamma distributions with shape parameters \(\alpha_k\) and scale parameters \(\beta\).

b) Provide some evidence that your method works by constructing some suitable sanity checks.

a) Construct an efficient algorithm for simulating realizations of a bivariate standard normal with zero correlation truncated to a disc with radius \(\rho\) centered at \(\delta,0\). Implement the algorithm in the form of a R function taking \(n, \rho, \delta\) as arguments that returns \(n\) realizations in \((n \times 2)\)-matrix. Test your function using a suitable sanity check. If using some form of rejection sampling, examine the overall acceptance rate and its dependency on \(\rho\) and \(\delta\). If needed you're allowed to use built-in quantile functions in R in addition to runif.

b) Building on the function implemented in point a), using a suitable linear transformation of the coordinates, construct a more general function simulating from a bivariate normal with mean vector \((\mu_x,\mu_y)\), zero correlation, both standard deviations equal to \(\sigma\), truncated to a disc with radius \(r\) centered at \(x_0,y_0\). Again, test your function using a suitable sanity check.

Let \(X=\frac{V_1/d_1}{V_2/d_2}\) where \(V_1\) and \(V_2\) are independent chi-square distributed random variables with \(d_1\) and \(d_2\) degrees of freedom, respectively.

a) Show that the density of \(X\) is \[ f_X(x)=\frac{\sqrt{\frac{(d_1 x)^{d_1} d_2^{d_2}}{(d_1 x+d_2)^{d_1+d_2}}}}{x\,\mathrm{B}\!\left(\frac{d_1}{2},\frac{d_2}{2}\right)} \] where \(B\) is the Beta function.

Hint: First, using the transformation formula for joint densities, derive the joint density of \(X\) and some strategically choosen function \(Y\) of \(V_1\) and \(V_2\). Next, integrate out \(y\) to find the (marginal) density of \(X\).

b) Verify that the result is correct by simulation. What is the name of this distribution?

a) Derive a method for simulating from the random variable \(Y\) having density (omitting the normalising constant) \[ f^*(y) = \begin{cases} \big(\sqrt{y^2 + 4} - y\big)^2 & \text{for }y\ge 0 \\ 0 & \text{for }y < 0. \end{cases} \]

Hint: Are \(f^*(y)\) and \(y^2 f^*(y)\) bounded? What does the limiting value of \(y^2 f^*(y)\) as \(y \rightarrow \infty \) tell you about the tail behaviour of the distribution?

b) Implement the method as a function in R and test that the method works. If your algorithm uses a form of rejection sampling, estimate or derive the overall acceptance probability.

c) According the strong law of large number, under certain conditions, \(\bar Y_n= \frac1n\sum_{i=1}^n Y_i\) where \(Y_1,Y_2,\dots\) are iid observations from \(f(y)\) should converge almost surely to \(E(Y)\) as \(n\rightarrow\infty\). Examine if the law appears to hold in this case by plotting a simulated sequence of \(\bar Y_n\). Comment on what you see.

Hint: You may want to use the function

cummean <- function(x) 
  cumsum(x)/(1:length(x))

to compute the sequence \(\bar Y_1,\bar Y_2,\dots\)

If \(X \sim N(\mu,\sigma^2)\) and \(Y = \frac 1{1+e^{-X}}\), then \(Y\) has a so-called logit-normal distribution with parameters \(\mu\) and \(\sigma\). This distribution appears in generalized linear mixed models with a binomial response, a logit link and the linear predictor contains normally distributed random effects and its expected value is of particular interest. A formula for the expected value is not available in closed form, however.

In this problem the aim is to estimate \(E(Y)\) using simulations instead. It clearly follows that \(E(Y)\) can be expressed as \[ E(h(X)) \tag{1} \] where \(X\sim N(\mu,\sigma^2)\) and \(h(X)=\frac 1{1+e^{-X}}\).

We will assume that \(\mu = 1\) and that \(\sigma = 3\).

a) First compute a Monte-Carlo estimate of (1) based on \(n=100\) realisations \(x_1,x_2,\dots,x_n\) of \(X\). Also compute an estimate of the standard error of the Monte-Carlo estimator.

b) Using the samples from \(f(x)\) in point a), construct \(n\) additional antithetic samples \(x_1^*,x_2^*,\dots,x_n^*\) from \(f(x)\) and compute a Monte-Carlo estimate of (1) using the combined \(2n\) samples. Explain how a standard error of this estimate of (1) can be computed. How does the standard error compare to the estimator in a)? Also estimate the correlation between \(h(X)\) and \(h(X^*)\) from the simulated realizations.

c) Instead of antithetic sampling, consider importance sampling. First make a graph of \(f(x)\) (the pdf of \(X\)), \(h(x)\), and the product between these two functions. Use a proposal density \(g(x)\) suggested by the graph to compute an improved importance sampling estimate of (1) based on \(n\) samples generated from \(g(x)\). One way of choosing the parameters of \(g(x)\) is perhaps to match the curvature of the mode to that of \(h(x)f(x)\). How does standard error of the importance sampling estimate compare to that of the estimators in point a) and b)?

d) Attempt to combine antithetic sampling with importance sampling and compare the standard error of the estimate with the other approaches. For your choice of \(g(x)\) in point c), make a graph of \(h(x)f(x)/g(x)\). Is it a monotonic function and why is that desirable?

e) Compare the above estimates to what you obtain by evaluating \[ \int_{-\infty}^\infty h(x) f(x) dx \] numerically using the integrate function in R.

Edit: Points b to e were modified on January 18

a) Write a function in R that simulates \(n\) realisations of \(\mathbf X\) from a multivariate normal distribution with mean vector \(\boldsymbol\mu\) and variance matrix \(\boldsymbol\Sigma\), either using the eigen decomposition or the Cholesky factor decomposition of the variance matrix.

b) A Gaussian copula is a continuous distribution over the unit square defined by \(U_1=\Phi(X_1)\) and \(U_2=\Phi(X_2)\) where \(\Phi\) is the standard normal cdf and \(X_1,X_2\) are jointly bivariate standard normal with correlation \(\rho\). Building on the function in point a), write an R function that simulates \(n\) realizations of \(U_1,U_2\). Generate a sample of \(1000\) realizations. Create histograms of \(U_1\) and \(U_2\). Are you able to prove the pattern that is indicated by the histogram? Also make a scatter plot of the realizations. Denoting the correlation between \(U_1,U_2\) by \(\gamma\), compute an estimate of \(\gamma\) when \(\rho=0.5\) using Monte-Carlo integration.

c) Are you able to derive an analytic relation between \(\gamma\) and \(\rho\)? Hints: This involves finding \(E(\Phi(X_1)\Phi(X_2))\). Note that \(\Phi(X_1)=E(I(Z_1\le X_1)|X_1)\) where \(Z_1\) is an independent standard normal random variable. Using the law of total expectation, prove that \(E(\Phi(X_1)\Phi(X_2))=P(Z_1\le X_1 \cap Z_2 \le X_2)\) where \(Z_2\) is another independent standard normal variable. This probability can in turn be re-expressed as the probability that independent standard normal random variables falls in a certain region in the plane. A closed form expression for this probability can then in turn be easily derived after a change to polar coordinates. Does your analytic formula match the simulation results in point b)?

d) Let \(Q(u;n,p)\) be the quantile function of a binomial distrution with parameters \(n\) and \(p\) (the function qbinom in R) and let \(Y_1 = Q(U_1;10, 0.25)\) and \(Y_2=Q(U_2;10, 0.75)\). Building on the R-function in point b) write an R-function that simulates \(n\) realisations of \(Y_1,Y_2\). What are the marginal distributions of \(Y_1\) and \(Y_2\)? Make a plot showing the simulated empirical joint distribution of \(n=100000\) realisations. Are \(Y_1\) and \(Y_2\) independent? Estimate their correlation by Monte-Carlo integration.