Computational-Statistics-732A90 Computer Lab 2 Solution

Linko¨pings Universitet, IDA, Statistik

Kurskod och namn: 732A90 Computational Statistics
Delmomentsansvarig: Krzysztof Bartoszek, Bayu Brahmantio, Jaskirat Marar, Shashi Nagarajan
Instruktioner: This computer laboratory is part of the examination for the
Computational Statistics course
Create a group report, (that is directly presentable, if you are a presenting group), on the solutions to the lab as a .PDF file.
Be concise and do not include unnecessary printouts and figures produced by the software and not required in the assignments.
All R code should be included as an appendix into your report. A typical lab report should 2-4 pages of text plus some amount of figures plus appendix with codes.
In the report reference ALL consulted sources and disclose ALL collaborations.
The report should be handed in via LISAM
The report has to be written in English.
Question 1: Optimizing parameters
Finding the minimum or maximum of a function is usually presented as a goal in itself. Here you are asked to use the function optim() to create a procedure to approximate another function, through so–called parabolic interpolation. For this exercise let f(x) be a continuous function on the interval [0,1] and let x0,x1,x2 ∈ [0,1] such that f(x1) < f(x0),f(x2). We will approximate the function f(x) with a function that is piecewise f˜(x) = a0 + a1x + a2x2, i.e. a piecewise quadratic function.
1. Write a function that uses optim() and finds values of (a0,a1,a2) for which f˜ interpolates f at user provided points x0,x1,x2. Interpolate means f(x0) = f˜(x0), f(x1) = f˜(x1) and f(x2) = f˜(x2). optim() should minimize the squared error, i.e. find (a0,a1,a2) that make (f(x0) − f˜(x0))2 + (f(x1) − f˜(x1))2 + (f(x2) − f˜(x2))2 as small as possible.
2. Now construct a function that approximates a function defined on the interval [0,1]. Your function should take as a parameter the number of equal–sized intervals that [0,1] is to be divided into and the function to approximate. The target function is known at the ends of the interval and also at the mid–point of the interval. Independently on each subinterval you should approximate the target function using the parabolic interpolater implemented in the previous part i.e. use the parabolic interpolater to find a0,a1,a2 for each subinterval.
3. Apply your function from the previous item to f1(x) = −x(1−x) and f2(x) = −xsin(10πx). Plot f1(·), f˜1(·) and f2(·), f˜2(·). How did your piecewise–parabolic interpolater fare? Explain what you observe. Take the number of subintervals to be at least 100.
Question 2: Maximizing likelihood
The file data.RData contains a sample from normal distribution with some parameters µ, σ. For this question read ?optim in detail.
1. Load the data to R environment.
2. Write down the log-likelihood function for 100 observations and derive maximum likelihoodestimators for µ, σ analytically by setting partial derivatives to zero. Use the derived formulae to obtain parameter estimates for the loaded data.
3. Optimize the minus log–likelihood function with initial parameters µ = 0, σ = 1. Try both Conjugate Gradient method (described in the presentation handout) and BFGS (discussed in the lecture) algorithm with gradient specified and without. Why it is a bad idea to maximize likelihood rather than maximizing log–likelihood?
4. Did the algorithms converge in all cases? What were the optimal values of parametersand how many function and gradient evaluations were required for algorithms to converge?
Which settings would you recommend?