STATS 782 Statistical Computing Assignment 4 (2023; Semester 1)
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
Department of Statistics
STATS 782 Statistical Computing
Assignment 4 (2023; Semester 1)
Total: 80 marks
Due: 10:00 NZDT, May 29, 2023
Notes:
(i) Please see Canvas Piazza for any additional instructions.
(ii) Remember to comment on all your output using your own words!
1. [20 marks] Use S3 to answer this question. We are interested in vectors on the plane. A
vector has a starting position and an ending position, where the former is the origin by default. Each vector must satisfy the following conditions:
❼ The x and y coordinates of both starting and end positions contains no NA, NaN or Inf values.
❼ If myvec1 contains n vectors then they are stored in a n × 4 matrix.
(a) Write a constructor function to create vectors. It can be a little like the "coords" example in the notes. So something like
n <- 8
set.seed(UPI.last) # UPI. last = last digit of your UPI
myvecs1 <- Vectors3(0 , 0 , rnorm(n), rnorm(n))
myvecs2 <- Vectors3(rnorm(n), rnorm(n), rnorm(n), rnorm(n))
should work. Use recycling when natural. The first produces 8 vectors emanating from the origin and potentially pointing in any direction. The second are 8 vectors starting ‘near’ the origin and potentially pointing in any direction too. [5 marks]
(b) Write a method to return the (Euclidean) length of the vectors. Then apply your function to myvecs1 and myvecs2 separately. [2 marks]
(c) Write a method to plot vectors. The user should be able to change some major characteris- tics such as the line width and colour. Then apply your function to myvecs1 and myvecs2 separately—show several different options to demonstrate that your function works quite flexibly. [9 marks]
(d) By increasing the value of n, estimate the asymptotic mean length of myvecs1 and myvecs2. Call them estimates of E(L1 ) and E(L2 ) respectively. For d = 2, can you
use simple results from mathematical statistics relating toE(L1(d)) and E(L2(d))—and demon-
strate it with your data? [4 marks]
2. [40 marks] In this question we we will look at polynomials, and later, a special type called
orthogonal polynomials (albeit, simplifed a lot). We will do this question using S4.
(a) Define polynomial objects using the class "poly4". It should carry information on its coefficients, and the domain upon which it is defined—assumed to be a simple interval on . of the form (a, b). Allow for some internal checking. Then write a constructor function; make the default the constant 1 defined over the whole real line. Hint: rather than storing the degree of the polynomial, delete any leading 0 coefficients so that its degree is simply the length of its coefficients vector minus one, e.g., 2x + 0x4 is really of degree 1. [4 marks]
(b) Use your code to construct objects representing the following polynomials over [ -1, 1]: P0 (x) = 1, P1 (x) = x, P2 
(c) Write three accessor functions: to extract the coefficients of your polynomials, to extract the domain of the polynomial, and return the degree of the polynomial. Run all three on P3 (x). [3 marks]
(d) Write a show() method to print out the individual objects, preferably nicely. However, just printing out the coefficients with a one sentence description explaining it will give you some marks. The following commands should work and print out something informative. Run your function on P0 (x) and P3 (x) like something as follows. [6 marks]
P0
P3
(e) Given polynomials
A2 (x) = a0 + a1 x + a2 x2 , (1)
B3 (x) = b0 + b1 x + b2 x2 + b3 x3 ,
C2 (x) = c0 + c1 x + c2 x2 ,
D3 (x) = d0 + d1 x + d2 x2 + d3 x3 ,
and a constant k, prove (vii), where the coefficients are written like a vector in ascending
order of degree. [3 marks]
(i) kA2 (x) = (ka0 , ka1 , ka2 ),
(ii) A2 (x) + B3 (x) = (a0 + b0 , a1 + b1 , a2 + b2 , b3 ),
(iii) A2 (x) + C2 (x) = (a0 + c0 , a1 + c1 , a2 + c2 ),
(iv) A2 (x) . B3 (x) = (a0 b0 , a1 b0 + a0 b1 , a0 b2 + a1 b1 + a2 b0 , a0 b3 + a1 b2 + a2 b1 , a1 b3 + a2 b2 , a2 b3 ), (v) A2 (x) . C2 (x) = (a0 c0 , a1 c0 + a0 c1 , a0 c2 + a1 c1 + a2 c0 , a1 c2 + a2 c1 , a2 c2 ),
(vi) B3 (x) . C2 (x) = (b0 c0 , b1 c0 + b0 c1 , b0 c2 + b1 c1 + b2 c0 , b1 c2 + b2 c1 + b3 c0 , b2 c2 ),
(vii) B3 (x) . D3 (x) = (b0 d0 , b1 d0 + b0 d1 , b0 d2 + b1 d1 + b2 d0 , b3 d0 + b2 d1 + b1 d2 + b0 d3 , b3 d1 + b2 d2 + b1 d3 , b3 d2 + b2 d3 , b3 d3 ).
Note: these formulas will be used below to obtain the pattern for selecting a subset of each polynomial’s coefficients and taking the inner product. The subscripts add up to the degree of the xd term. If the coefficients are sorted in ascending order then the computation is not difficult.
(f) Ideally we want to allow for ‘addition’ and ‘multiplication’. Let’s do it first for the simple case of ‘addition’ (and ‘subtraction’ of course) only, where one of the operands is a single numeric. Implement this and then run it (your names for P0 (x) etc. will be different) on
2 + P2
-1 - P2
P3 + 1
P3 - 3
1:2 + P3 # This should give an error
That is, write some S4 code to allow addition (and subtraction). Do this for case of: ((1)) one polynomial and a scalar numeric (it should give an error if a vector of length > 1 is used). [5 marks]
(g) Now we want to allow for ‘addition’ (and ‘subtraction’) between two "poly4" objects, e.g.,
P0 + P0
P1 + P1
P2 + P3
P3 - P2
Implement this and then run it. Hint: use the formulas based on (1). [6 marks]
(h) Now write some S4 code to allow multiplication of your polynomials. We need to do this for two cases: between ((1)) one polynomial and a scalar numeric (it should give an error if a vector of length > 1 is used); and ((2)) two polynomials. Let’s do the first case first. [4 marks]
2 * P3
P2 * 10
(1:2) * P3 # This should give an error
(i) Now implement multiplication between two polynomials. For this, you’ll need to check that both polynomials have the same domain. Run your code on
P0 * P0
P1 * P1
P1 * P2
P2 * P3
Hint: as it can be done without a for() loop, more marks may be awarded if you achieve this. [7 marks]
3. [20 marks] This question is a continuation of Q2.
Now we are interested in a special type of polynomial called an orthogonal polynomial. We say that Pm (x) and Pn (x) are orthogonal with respect to the (weight) function w(x) over the interval (a, b) if
la b Pm (x) Pn (x) w(x) dx = 0. (2)
(a) Use setClass() to create a subclass "opoly4" of the previously defined S4 class "poly4" that contains information about w(x) too. Don’t worry about a validity function yet, however, make the 0th degree Legendre polynomial P0 on [-1, 1], the default. [4 marks]
(b) The most well-known example of a sequence of orthogonal polynomials are called Legendre polynomials, which have w(x) = 1 defined on [-1, 1]. The first few Legendre polynomials are P0 (x), P1 (x), P2 (x), P3 (x) given earlier, and defined on [-1, 1]. Write a constructor function and then create objects of these using your code. [3 marks]
(c) Write a generic function integ(object1, object2) which returns the quantity (2) over its domain. Test your code like the following.
Dim <- 3
matinteg <- matrix(NA_real_ , Dim + 1, Dim + 1)
alllist <- list(P0, P1, P2, P3) # A sequence of orthog polys
for (i in 0:Dim)
for (j in i:Dim)
matinteg[i + 1, j + 1] <-
matinteg[j + 1, i + 1] <- integ(alllist[[i + 1]],
alllist[[j + 1]])
round(matinteg, 5)
Hint: do some internal error checking, and build your code around integrate(). [10 marks]
(d) Write a generic (why?) function is. orthog(object1, object2) which tests whether the two polynomials are orthogonal. Test it at least on a few examples of pairs of orthogonal polynomials and nonorthogonal polynomials, e.g.,
is.orthog(newopoly4(1:4),
newopoly4(4:1)) # Should be FALSE
is.orthog(newopoly4(c(0 , 0 , 1)),
newopoly4(1:3)) # Should be FALSE
# Try this test for Legendre polynomials:
matorthog <- matrix(FALSE, Dim + 1, Dim + 1)
for (i in 0:Dim)
for (j in i:Dim)
matorthog[i + 1, j + 1] <-
matorthog[j + 1, i + 1] <- is.orthog(alllist[[i + 1]],
alllist[[j + 1]])
matorthog
Comment: is. orthog() could be incorporated into a validity function for "opoly4" objects, but we won’t do that—because actually the Legendre polynomials are a sequence that satisfies a recurrence formula—but this is too complicated for 782. [3 marks]
2023-06-20