Math 365 / Comp 365: Final Exam
Monday, May 4 - Monday, May 11, 2015
Upload an html file to Moodle and turn in any hand written work by noon on Monday, May 11, 2015
Tyler J. Skluzacek
I pledge that all work on this exam is my own, and that I have not consulted with other people.
The following line sources functions from the class file 365Functions.r. Feel free to use any of these functions.
source("http://www.macalester.edu/~dshuman1/data/365/365Functions.r")
source("http://www.macalester.edu/~dshuman1/data/365/2015.365.Final.Data.r")
require(Matrix)
require(expm)Problem 1
Part A
Here is an image of the transition matrix. I’ve attached a picture of the Markov Chain on an additional document.
CnL = spMatrix(20,20, i = seq(1,19), j = seq(2,20), x=rep(1/3,19))+spMatrix(20,20, i = seq(1,18), j=seq(3,20), rep(1/3,18))+spMatrix(20,20,i=seq(1,20),j=seq(1,20),x=rep(1/3,20))
CnL[,6] = 0
CnL[,8] = 0
CnL[,9] = 0
CnL[,18] = 0
CnL[,20] = 0
CnL[6,] = 0
CnL[8,] = 0
CnL[9,] = 0
CnL[18,] = 0
CnL[20,] = 0
CnL[4,15] = 1/3
CnL[5,15] = 1/3
CnL[7,17] = 1/3
CnL[7,3] = 1/3
CnL[16,12] = 1/3
CnL[17,12] = 1/3
CnL[19,20] = 2/3
CnL[20,20] = 1
image(CnL)
CnL = as.matrix(CnL) #To save entry-edits. Part B
tCNL = as.matrix(t(CnL))
x20 = cbind(c(1, rep(0,19)))
turn20 = tCNL%^%20%*%x20
turn20[20,]*100## [1] 57.62Therefore, the percentage chance that I’m on square 20 after 20 turns is 57.6%.
Part C
cnl2 = CnL
cnl2[20,20] = 1/3
cnl2[20,1] = 1/3
cnl2[20,2] = 1/3
cnl2[19,1] = 1/3
cnl2[19,20] = 1/3
image(cnl2) #New looping matrix. 
tCNL2 = t(as.matrix(cnl2))
turnProb = tCNL2%^%1000000%*%x20
turnProb## [,1]
## [1,] 0.04850
## [2,] 0.04042
## [3,] 0.05783
## [4,] 0.04912
## [5,] 0.05347
## [6,] 0.00000
## [7,] 0.02674
## [8,] 0.00000
## [9,] 0.00000
## [10,] 0.00000
## [11,] 0.00000
## [12,] 0.11808
## [13,] 0.05904
## [14,] 0.08856
## [15,] 0.12510
## [16,] 0.10683
## [17,] 0.12933
## [18,] 0.00000
## [19,] 0.06467
## [20,] 0.03233max(turnProb) ## [1] 0.1293After 1 million turns, it is safe to assume that the probabilities of being on a given space has converged (because one would probably stop playing by 1,000,000 turns anyways).
One is most likely to be on square 17 (~12.9% probability).
The probability of being on square 10 is 0% (because 8 and 9 are unreachable nodes).
The probability of being on square 19 is ~6.4%.
Problem 2
Wolfram Alpha was used extensively in integral calculation.
Part A.
To begin solving this problem, I first need to write out the polynomial equation using the degree-2 Legendre Polynomial:
\[f(x) = \frac{1}{\sqrt(2)}a_0 + \sqrt(\frac{3}{2})xa_1 + \sqrt(\frac{5}{8})(3x^2 - 1)a_2\]
Now I can find the inner product between h (\(e^x\)) and each of the Legendre terms in order to find the terms of the approximating polynomial.
\[a_0 = <h,q_0> = \int_{-1}^1 \sqrt(\frac{1}{2})e^x dx\]
\[a_1 = <h,q_1> = \int_{-1}^1 \sqrt(\frac{3}{2})xe^x dx\]
\[a_2 = <h,q_2> = \int_{-1}^1 \sqrt(\frac{5}{8})(3x^2 - 1)e^x dx\]
Simplifying this, I get:
\[a_0 = 1.662\] \[a_1 = 0.901\] \[a_2 = 0.226\]
Therefore, a degree-2 approximation can be written as the following:
\[f(x) = \sqrt(\frac{1}{2})(1.662) + \sqrt(\frac{3}{2})x(0.901) + \sqrt(\frac{5}{8})(3x^2 - 1)(0.226)\].
Part B.
First, I can write the equation for the functional approximation’s coefficients as follows (working backwards from the previous method in Part A):
\[a_0q_0+a_1q_1+a_2q_2=2x^2-4x+7\]
Now, for each term in the above sequence, I can set up the approximation using the Legendre Terms derived in Part A as follows, such that g(x) minimizes \(||g-f||\).
\[\frac{1}{\sqrt(2)}\int_{-1}^1g(x)\frac{1}{\sqrt(2)}dx+\sqrt(\frac{3}{2})x\int_{-1}^1g(x)\sqrt(\frac{3}{2})xdx+\sqrt(\frac{5}{8})(3x^2-1)\int_{-1}^1g(x)\sqrt(\frac{5}{8})(3x^2-1)dx=2x^2-4x+7\]
Now substituting the Legendre Polynomials for a degree-2 approximation, we can remove the coefficients on the \(x\) and \(g(x)\) terms, and bring them outside the integral for simpler solving later.
\[\frac{1}{2}\int_{-1}^1g(x)dx+\frac{3}{2}x\int_{-1}^1g(x)xdx+\frac{5}{8}(3x^2-1)\int_{-1}^1g(x)(3x^2-1)dx=2x^2-4x+7\]
Now, we have a formula for all degree-3 polynomials with a degree-2 approximation of \(2x^2-4x+7\). One just needs to solve for \(g(x)\).
Problem 3
Here is the digram frequency matrix of Lincoln’s Gettysburg address:
print(G)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 0 1 2 5 0 1 4 0 2 0 1 9 0
## [2,] 1 0 0 0 5 0 0 0 1 0 0 1 0
## [3,] 12 0 0 0 4 0 0 2 1 0 0 0 0
## [4,] 2 0 1 6 14 0 0 3 13 0 0 0 0
## [5,] 16 3 8 26 3 5 2 7 6 0 0 4 5
## [6,] 3 0 0 1 0 1 0 0 5 0 0 0 0
## [7,] 5 1 0 0 5 0 1 4 0 0 0 1 0
## [8,] 24 0 0 0 32 1 0 0 7 0 0 1 0
## [9,] 0 1 8 1 3 0 2 0 0 0 0 2 0
## [10,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [11,] 0 0 0 0 1 0 0 0 0 0 0 0 0
## [12,] 3 0 0 4 6 0 0 1 6 0 0 8 2
## [13,] 2 1 0 0 7 0 0 0 1 0 0 0 0
## [14,] 10 0 5 9 4 1 9 0 2 0 0 3 2
## [15,] 1 3 1 3 0 6 1 1 0 0 0 1 4
## [16,] 0 0 0 0 5 0 0 0 0 0 0 4 0
## [17,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [18,] 8 0 0 1 26 4 3 0 1 0 1 3 0
## [19,] 4 2 2 0 10 2 1 6 1 0 1 0 0
## [20,] 4 1 4 1 11 5 1 47 18 0 0 3 0
## [21,] 1 0 0 0 0 0 3 0 0 0 0 2 0
## [22,] 2 0 0 0 17 0 0 0 3 0 0 0 0
## [23,] 2 1 0 0 11 0 0 8 1 0 0 0 0
## [24,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [25,] 2 0 0 1 1 0 1 1 0 0 0 0 0
## [26,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 15 0 1 0 10 5 36 1 8 0 0
## [2,] 0 1 0 0 2 0 0 2 0 0 0
## [3,] 0 7 0 0 4 0 1 0 0 0 0
## [4,] 0 4 1 0 0 4 4 1 1 4 0
## [5,] 10 5 4 1 22 9 12 2 4 8 0
## [6,] 0 10 0 0 3 0 3 1 0 0 0
## [7,] 0 3 1 0 6 0 0 0 0 1 0
## [8,] 0 8 0 0 0 0 5 1 0 0 0
## [9,] 16 9 0 0 2 9 8 0 7 0 0
## [10,] 0 0 0 0 0 0 0 0 0 0 0
## [11,] 0 0 0 0 1 0 0 0 0 1 0
## [12,] 3 3 0 0 1 0 1 0 1 1 0
## [13,] 0 0 0 0 0 0 2 0 0 0 0
## [14,] 4 12 0 0 0 4 8 1 1 0 0
## [15,] 20 2 5 0 17 3 13 7 2 3 0
## [16,] 0 4 0 0 2 0 0 0 0 0 0
## [17,] 0 0 0 0 0 0 0 1 0 0 0
## [18,] 1 6 2 0 0 5 12 3 0 3 0
## [19,] 1 4 0 0 1 0 8 1 0 0 0
## [20,] 2 11 1 0 2 0 9 0 0 5 0
## [21,] 3 0 0 0 5 5 2 0 0 0 0
## [22,] 0 2 0 0 0 0 0 0 0 0 0
## [23,] 1 2 0 0 0 0 1 0 0 1 0
## [24,] 0 0 0 0 0 0 0 0 0 0 0
## [25,] 1 0 0 0 1 0 1 0 0 1 0
## [26,] 0 0 0 0 0 0 0 0 0 0 0
## [,25] [,26]
## [1,] 1 0
## [2,] 1 0
## [3,] 0 0
## [4,] 0 0
## [5,] 3 0
## [6,] 0 0
## [7,] 0 0
## [8,] 0 0
## [9,] 0 0
## [10,] 0 0
## [11,] 0 0
## [12,] 2 0
## [13,] 0 0
## [14,] 2 0
## [15,] 0 0
## [16,] 0 0
## [17,] 0 0
## [18,] 0 0
## [19,] 0 0
## [20,] 1 0
## [21,] 0 0
## [22,] 0 0
## [23,] 0 0
## [24,] 0 0
## [25,] 0 0
## [26,] 0 0Here is the vector \(f=A1\) for part (a):
letters = c('A','B','C','D','E','F','G','H','I','J','K','L','M',
'N','O','P','Q','R','S','T','U','V','W','X','Y','Z')
f = G%*%rep(1,26)
t = cbind(letters,f)
print(t)## letters
## [1,] "A" "102"
## [2,] "B" "14"
## [3,] "C" "31"
## [4,] "D" "58"
## [5,] "E" "165"
## [6,] "F" "27"
## [7,] "G" "28"
## [8,] "H" "79"
## [9,] "I" "68"
## [10,] "J" "0"
## [11,] "K" "3"
## [12,] "L" "42"
## [13,] "M" "13"
## [14,] "N" "77"
## [15,] "O" "93"
## [16,] "P" "15"
## [17,] "Q" "1"
## [18,] "R" "79"
## [19,] "S" "44"
## [20,] "T" "126"
## [21,] "U" "21"
## [22,] "V" "24"
## [23,] "W" "28"
## [24,] "X" "0"
## [25,] "Y" "10"
## [26,] "Z" "0"Here is the scrambled matrix for part (e):
print(S)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 0 3 0 4 0 0 0 4 0 3 2 0 0
## [2,] 6 0 24 23 0 7 0 2 21 1 17 3 2
## [3,] 0 30 2 20 0 1 0 13 4 32 11 3 1
## [4,] 2 12 7 0 0 10 5 0 0 33 3 8 0
## [5,] 0 0 0 0 0 0 0 0 0 0 0 0 0
## [6,] 0 8 0 10 0 12 0 1 0 8 0 3 0
## [7,] 0 0 0 1 0 0 0 0 0 0 0 0 0
## [8,] 0 10 0 47 0 0 0 1 2 10 1 29 0
## [9,] 0 5 2 7 0 1 0 8 1 27 3 4 0
## [10,] 3 9 19 8 0 15 0 26 4 13 14 3 3
## [11,] 0 14 0 2 0 0 0 0 7 5 0 0 0
## [12,] 0 70 0 52 0 2 0 7 2 24 0 2 0
## [13,] 0 2 0 0 0 0 0 0 0 3 0 0 0
## [14,] 0 19 4 18 0 4 0 2 4 25 4 4 4
## [15,] 0 15 2 18 0 1 0 5 0 17 0 20 0
## [16,] 0 1 0 1 0 0 0 0 0 0 0 0 0
## [17,] 3 36 5 36 0 11 0 15 3 62 4 21 0
## [18,] 7 27 7 53 1 2 0 3 16 59 3 4 0
## [19,] 17 30 0 13 0 0 0 1 3 29 0 21 0
## [20,] 0 27 1 21 0 3 0 1 5 13 2 0 0
## [21,] 0 32 3 88 0 2 0 12 22 102 7 202 0
## [22,] 0 8 2 7 0 7 0 0 0 0 3 4 0
## [23,] 4 34 79 21 2 1 0 14 4 46 3 3 4
## [24,] 3 36 10 42 0 2 0 2 7 26 19 3 2
## [25,] 0 7 0 8 0 0 0 0 1 4 0 0 0
## [26,] 1 86 80 25 3 13 0 33 7 27 8 20 2
## [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 0 2 0 2 1 0 2 12 0 2 1
## [2,] 1 14 0 35 29 25 33 93 5 102 53
## [3,] 3 5 0 26 26 2 1 19 7 2 4
## [4,] 11 24 1 77 15 43 10 64 0 111 55
## [5,] 0 0 0 0 0 0 0 0 6 0 0
## [6,] 1 0 0 3 15 0 0 4 2 0 12
## [7,] 0 0 0 0 0 0 0 0 0 0 0
## [8,] 0 0 0 1 2 2 2 3 0 4 0
## [9,] 3 0 0 8 1 3 8 18 0 3 2
## [10,] 100 10 0 19 48 10 62 50 52 74 21
## [11,] 0 0 0 0 8 0 0 0 8 0 17
## [12,] 1 0 0 3 2 0 0 24 4 0 2
## [13,] 0 0 0 0 0 0 0 0 11 0 0
## [14,] 12 8 0 3 48 7 4 24 8 4 4
## [15,] 12 3 0 1 25 2 1 6 5 5 5
## [16,] 0 0 0 0 0 0 0 2 1 0 0
## [17,] 14 2 0 28 1 13 4 73 11 11 16
## [18,] 5 8 0 21 4 10 9 20 11 15 4
## [19,] 1 0 0 2 9 0 1 4 5 0 3
## [20,] 0 0 0 5 11 0 2 19 10 4 3
## [21,] 22 1 0 14 23 1 3 21 14 6 9
## [22,] 4 4 1 36 26 1 3 18 0 27 13
## [23,] 1 53 0 31 1 11 2 63 4 12 6
## [24,] 7 2 0 12 2 2 4 16 2 5 63
## [25,] 0 0 0 0 0 0 0 0 0 0 0
## [26,] 18 30 3 70 100 36 28 89 4 71 26
## [,25] [,26]
## [1,] 0 8
## [2,] 25 0
## [3,] 0 35
## [4,] 11 23
## [5,] 0 0
## [6,] 0 15
## [7,] 0 4
## [8,] 0 36
## [9,] 2 7
## [10,] 4 2
## [11,] 1 42
## [12,] 0 162
## [13,] 0 2
## [14,] 0 6
## [15,] 0 23
## [16,] 0 0
## [17,] 0 28
## [18,] 2 106
## [19,] 0 29
## [20,] 0 52
## [21,] 2 56
## [22,] 0 6
## [23,] 1 58
## [24,] 2 50
## [25,] 0 56
## [26,] 26 46Part A
f tells us about the number of each letter that appears in the text. I confirmed this by downloading the textfile and using MSword’s Find function.
Part B
The numerical rank is calculated by counting all singular values greater than \(E_machine * (largest singular value)\):
svdG = svd(G)
tol=(max(svdG$d))*.Machine$double.eps*norm(G,type="I")
sings=svdG$d
numerical.rank = sum(sings>tol)
numerical.rank## [1] 23So the numerical rank of this matrix is 23. This makes sense given the data because there are only 23 letters that occur in the Gettysburg Address. The numerical rank is the rank given machine error.
Part C
The best rank-1 approximation is given by zeroing out the r-k trailing singular values of G (https://inst.eecs.berkeley.edu/~ee127a/book/login/l_svd_low_rank.html$})
newSigma = spMatrix(26,26,i=seq(1,26),j=seq(1,26), x = rep(0,26))
newSigma[1,1] = max(svdG$d)
newSigma = as.matrix(newSigma)
image(newSigma)
G1 = (svdG$u)%*%newSigma%*%t(svdG$v)
G1 #Best Rank 1 approximation## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 8.62376 1.182799 3.043517 6.05736 12.09458 2.739895 2.143502
## [2,] 1.23758 0.169741 0.436768 0.86928 1.73567 0.393196 0.307609
## [3,] 3.15764 0.433089 1.114403 2.21794 4.42851 1.003229 0.784857
## [4,] 5.29126 0.725728 1.867405 3.71660 7.42085 1.681112 1.315184
## [5,] 11.56130 1.585700 4.080242 8.12070 16.21440 3.673195 2.873651
## [6,] 2.30097 0.315592 0.812064 1.61621 3.22705 0.731052 0.571924
## [7,] 2.54075 0.348479 0.896687 1.78463 3.56333 0.807233 0.631522
## [8,] 9.11997 1.250857 3.218641 6.40590 12.79050 2.897548 2.266839
## [9,] 4.81841 0.660874 1.700526 3.38447 6.75769 1.530881 1.197654
## [10,] 0.00000 0.000000 0.000000 0.00000 0.00000 0.000000 0.000000
## [11,] 0.25929 0.035564 0.091511 0.18213 0.36365 0.082381 0.064449
## [12,] 3.16500 0.434099 1.117001 2.22311 4.43883 1.005569 0.786687
## [13,] 1.59550 0.218832 0.563087 1.12068 2.23764 0.506913 0.396573
## [14,] 5.69595 0.781233 2.010228 4.00086 7.98841 1.809687 1.415772
## [15,] 6.27630 0.860832 2.215047 4.40850 8.80234 1.994073 1.560023
## [16,] 1.37444 0.188513 0.485072 0.96541 1.92762 0.436681 0.341628
## [17,] 0.01923 0.002637 0.006786 0.01351 0.02697 0.006109 0.004779
## [18,] 7.76682 1.065265 2.741084 5.45544 10.89275 2.467632 1.930503
## [19,] 4.42629 0.607092 1.562138 3.10904 6.20775 1.406299 1.100190
## [20,] 11.71447 1.606709 4.134299 8.22829 16.42922 3.721860 2.911723
## [21,] 1.40210 0.192306 0.494833 0.98484 1.96641 0.445468 0.348503
## [22,] 3.06886 0.420912 1.083068 2.15558 4.30399 0.975021 0.762788
## [23,] 3.18393 0.436695 1.123680 2.23640 4.46537 1.011582 0.791391
## [24,] 0.00000 0.000000 0.000000 0.00000 0.00000 0.000000 0.000000
## [25,] 0.89311 0.122496 0.315200 0.62733 1.25257 0.283755 0.221990
## [26,] 0.00000 0.000000 0.000000 0.00000 0.00000 0.000000 0.000000
## [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 9.05692 6.26758 4.328e-26 0.2596121 3.26527 1.255025 6.53211
## [2,] 1.29974 0.89945 6.211e-27 0.0372563 0.46859 0.180106 0.93741
## [3,] 3.31625 2.29491 1.585e-26 0.0950586 1.19560 0.459535 2.39177
## [4,] 5.55704 3.84559 2.656e-26 0.1592896 2.00347 0.770043 4.00789
## [5,] 12.14201 8.40253 5.802e-26 0.3480446 4.37753 1.682529 8.75717
## [6,] 2.41655 1.67230 1.155e-26 0.0692691 0.87123 0.334863 1.74288
## [7,] 2.66837 1.84657 1.275e-26 0.0764874 0.96202 0.369758 1.92450
## [8,] 9.57805 6.62822 4.577e-26 0.2745501 3.45315 1.327239 6.90797
## [9,] 5.06044 3.50193 2.418e-26 0.1450549 1.82443 0.701229 3.64973
## [10,] 0.00000 0.00000 0.000e+00 0.0000000 0.00000 0.000000 0.00000
## [11,] 0.27232 0.18845 1.301e-27 0.0078058 0.09818 0.037735 0.19640
## [12,] 3.32398 2.30027 1.588e-26 0.0952802 1.19839 0.460607 2.39735
## [13,] 1.67564 1.15958 8.008e-27 0.0480313 0.60411 0.232195 1.20852
## [14,] 5.98205 4.13971 2.859e-26 0.1714725 2.15670 0.828938 4.31443
## [15,] 6.59155 4.56150 3.150e-26 0.1889435 2.37644 0.913397 4.75402
## [16,] 1.44348 0.99892 6.898e-27 0.0413766 0.52041 0.200024 1.04108
## [17,] 0.02019 0.01397 9.650e-29 0.0005788 0.00728 0.002798 0.01456
## [18,] 8.15694 5.64477 3.898e-26 0.2338145 2.94080 1.130314 5.88302
## [19,] 4.64862 3.21695 2.221e-26 0.1332504 1.67596 0.644164 3.35272
## [20,] 12.30288 8.51385 5.879e-26 0.3526557 4.43553 1.704820 8.87319
## [21,] 1.47253 1.01902 7.037e-27 0.0422092 0.53089 0.204049 1.06203
## [22,] 3.22300 2.23039 1.540e-26 0.0923857 1.16198 0.446614 2.32452
## [23,] 3.34386 2.31402 1.598e-26 0.0958500 1.20555 0.463361 2.41169
## [24,] 0.00000 0.00000 0.000e+00 0.0000000 0.00000 0.000000 0.00000
## [25,] 0.93797 0.64910 4.482e-27 0.0268865 0.33817 0.129976 0.67649
## [26,] 0.00000 0.00000 0.000e+00 0.0000000 0.00000 0.000000 0.00000
## [,15] [,16] [,17] [,18] [,19] [,20] [,21]
## [1,] 6.86963 1.513142 0.1441837 6.68589 3.730876 11.71052 1.600181
## [2,] 0.98584 0.217148 0.0206915 0.95948 0.535410 1.68055 0.229638
## [3,] 2.51536 0.554046 0.0527937 2.44808 1.366084 4.28788 0.585916
## [4,] 4.21498 0.928415 0.0884665 4.10225 2.289146 7.18520 0.981820
## [5,] 9.20965 2.028569 0.1932975 8.96333 5.001738 15.69952 2.145256
## [6,] 1.83294 0.403733 0.0384708 1.78391 0.995464 3.12457 0.426957
## [7,] 2.02394 0.445805 0.0424797 1.96981 1.099198 3.45018 0.471448
## [8,] 7.26490 1.600208 0.1524800 7.07060 3.945550 12.38435 1.692255
## [9,] 3.83831 0.845449 0.0805608 3.73566 2.084579 6.54310 0.894080
## [10,] 0.00000 0.000000 0.0000000 0.00000 0.000000 0.00000 0.000000
## [11,] 0.20655 0.045496 0.0043352 0.20103 0.112178 0.35210 0.048113
## [12,] 2.52122 0.555338 0.0529169 2.45379 1.369269 4.29788 0.587282
## [13,] 1.27096 0.279949 0.0266757 1.23697 0.690256 2.16659 0.296052
## [14,] 4.53735 0.999423 0.0952326 4.41600 2.464225 7.73474 1.056912
## [15,] 4.99966 1.101252 0.1049357 4.86594 2.715301 8.52282 1.164599
## [16,] 1.09487 0.241162 0.0229798 1.06559 0.594622 1.86641 0.255035
## [17,] 0.01532 0.003374 0.0003215 0.01491 0.008318 0.02611 0.003568
## [18,] 6.18699 1.362782 0.1298562 6.02152 3.360141 10.54685 1.441171
## [19,] 3.52595 0.776646 0.0740048 3.43165 1.914937 6.01063 0.821321
## [20,] 9.33167 2.055445 0.1958584 9.08209 5.068004 15.90752 2.173678
## [21,] 1.11690 0.246015 0.0234422 1.08703 0.606587 1.90396 0.260167
## [22,] 2.44463 0.538468 0.0513093 2.37925 1.327672 4.16732 0.569442
## [23,] 2.53630 0.558659 0.0532333 2.46846 1.377457 4.32358 0.590794
## [24,] 0.00000 0.000000 0.0000000 0.00000 0.000000 0.00000 0.000000
## [25,] 0.71145 0.156707 0.0149323 0.69242 0.386385 1.21279 0.165722
## [26,] 0.00000 0.000000 0.0000000 0.00000 0.000000 0.00000 0.000000
## [,22] [,23] [,24] [,25] [,26]
## [1,] 2.190809 2.79854 0 0.922642 0
## [2,] 0.314398 0.40161 0 0.132406 0
## [3,] 0.802178 1.02470 0 0.337831 0
## [4,] 1.344210 1.71709 0 0.566104 0
## [5,] 2.937072 3.75181 0 1.236925 0
## [6,] 0.584547 0.74670 0 0.246177 0
## [7,] 0.645460 0.82451 0 0.271830 0
## [8,] 2.316868 2.95956 0 0.975731 0
## [9,] 1.224086 1.56365 0 0.515514 0
## [10,] 0.000000 0.00000 0 0.000000 0
## [11,] 0.065872 0.08414 0 0.027741 0
## [12,] 0.804049 1.02709 0 0.338619 0
## [13,] 0.405326 0.51776 0 0.170700 0
## [14,] 1.447019 1.84842 0 0.609400 0
## [15,] 1.594453 2.03675 0 0.671491 0
## [16,] 0.349168 0.44603 0 0.147049 0
## [17,] 0.004885 0.00624 0 0.002057 0
## [18,] 1.973109 2.52045 0 0.830959 0
## [19,] 1.124471 1.43640 0 0.473562 0
## [20,] 2.975984 3.80152 0 1.253312 0
## [21,] 0.356194 0.45500 0 0.150008 0
## [22,] 0.779623 0.99589 0 0.328332 0
## [23,] 0.808857 1.03323 0 0.340644 0
## [24,] 0.000000 0.00000 0 0.000000 0
## [25,] 0.226890 0.28983 0 0.095553 0
## [26,] 0.000000 0.00000 0 0.000000 0newSigma2 = newSigma
newSigma2[2,2] = sort(svdG$d, partial = 25)[25]
image(newSigma2)
G2 = (svdG$u)%*%newSigma2%*%t(svdG$v)
G2 #Best Rank 2 approximation## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 4.314e+00 2.234e+00 5.187e+00 1.194e+01 -2.389e+00 3.477e+00
## [2,] 1.588e+00 8.434e-02 2.627e-01 3.911e-01 2.912e+00 3.333e-01
## [3,] 3.769e+00 2.839e-01 8.103e-01 1.383e+00 6.483e+00 8.987e-01
## [4,] 6.510e+00 4.284e-01 1.261e+00 2.052e+00 1.152e+01 1.473e+00
## [5,] 8.763e+00 2.268e+00 5.472e+00 1.194e+01 6.812e+00 4.152e+00
## [6,] 2.296e+00 3.168e-01 8.146e-01 1.623e+00 3.210e+00 7.319e-01
## [7,] 2.983e+00 2.406e-01 6.768e-01 1.181e+00 5.049e+00 7.316e-01
## [8,] 1.241e+01 4.474e-01 1.581e+00 1.907e+00 2.386e+01 2.334e+00
## [9,] 3.069e+00 1.088e+00 2.570e+00 5.773e+00 8.804e-01 1.830e+00
## [10,] 1.401e-23 -3.417e-24 -6.967e-24 -1.913e-23 4.707e-23 -2.395e-24
## [11,] 2.878e-01 2.860e-02 7.732e-02 1.432e-01 4.595e-01 7.750e-02
## [12,] 3.419e+00 3.721e-01 9.906e-01 1.876e+00 5.293e+00 9.621e-01
## [13,] 2.134e+00 8.748e-02 2.953e-01 3.852e-01 4.047e+00 4.149e-01
## [14,] 5.150e+00 9.145e-01 2.282e+00 4.747e+00 6.153e+00 1.903e+00
## [15,] 2.978e+00 1.665e+00 3.855e+00 8.913e+00 -2.279e+00 2.558e+00
## [16,] 1.701e+00 1.087e-01 3.224e-01 5.187e-01 3.027e+00 3.808e-01
## [17,] 1.108e-02 4.625e-03 1.084e-02 2.464e-02 -4.222e-04 7.502e-03
## [18,] 9.202e+00 7.152e-01 2.027e+00 3.495e+00 1.572e+01 2.222e+00
## [19,] 5.104e+00 4.417e-01 1.225e+00 2.183e+00 8.487e+00 1.290e+00
## [20,] 1.491e+01 8.267e-01 2.544e+00 3.861e+00 2.718e+01 3.175e+00
## [21,] 6.751e-01 3.696e-01 8.564e-01 1.978e+00 -4.765e-01 5.697e-01
## [22,] 4.668e+00 3.097e-02 2.881e-01 -2.790e-02 9.676e+00 7.017e-01
## [23,] 4.505e+00 1.144e-01 4.667e-01 4.320e-01 8.905e+00 7.857e-01
## [24,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## [25,] 8.555e-01 1.317e-01 3.339e-01 6.787e-01 1.126e+00 2.902e-01
## [26,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 3.835e+00 1.057e-01 1.387e+00 4.328e-26 3.917e-01 5.582e+00
## [2,] 1.702e-01 2.027e+00 1.296e+00 6.211e-27 2.653e-02 2.804e-01
## [3,] 5.448e-01 4.586e+00 2.987e+00 1.585e-26 7.632e-02 8.669e-01
## [4,] 8.368e-01 8.088e+00 5.226e+00 2.656e-26 1.219e-01 1.348e+00
## [5,] 3.972e+00 6.331e+00 5.234e+00 5.802e-26 4.338e-01 5.882e+00
## [6,] 5.739e-01 2.406e+00 1.667e+00 1.155e-26 6.942e-02 8.739e-01
## [7,] 4.580e-01 3.587e+00 2.347e+00 1.275e-26 6.294e-02 7.244e-01
## [8,] 9.740e-01 1.642e+01 1.036e+01 4.577e-26 1.736e-01 1.683e+00
## [9,] 1.884e+00 1.428e+00 1.521e+00 2.418e-26 1.986e-01 2.765e+00
## [10,] -5.498e-24 2.909e-23 1.586e-23 0.000e+00 -4.292e-25 -7.530e-24
## [11,] 5.325e-02 3.316e-01 2.208e-01 1.301e-27 6.932e-03 8.284e-02
## [12,] 6.869e-01 3.852e+00 2.588e+00 1.588e-26 8.749e-02 1.062e+00
## [13,] 1.852e-01 2.794e+00 1.769e+00 8.008e-27 3.153e-02 3.147e-01
## [14,] 1.630e+00 4.847e+00 3.521e+00 2.859e-26 1.882e-01 2.450e+00
## [15,] 2.854e+00 -2.574e-01 8.273e-01 3.150e-26 2.900e-01 4.149e+00
## [16,] 2.133e-01 2.123e+00 1.369e+00 6.898e-27 3.136e-02 3.446e-01
## [17,] 7.978e-03 3.266e-03 4.745e-03 9.650e-29 8.285e-04 1.166e-02
## [18,] 1.367e+00 1.114e+01 7.270e+00 3.898e-26 1.898e-01 2.169e+00
## [19,] 8.340e-01 6.057e+00 3.985e+00 2.221e-26 1.125e-01 1.311e+00
## [20,] 1.657e+00 1.894e+01 1.213e+01 5.879e-26 2.547e-01 2.717e+00
## [21,] 6.338e-01 -3.730e-02 1.958e-01 7.037e-27 6.448e-02 9.217e-01
## [22,] 1.353e-01 6.543e+00 4.041e+00 1.540e-26 4.341e-02 3.027e-01
## [23,] 2.729e-01 6.088e+00 3.810e+00 1.598e-26 5.538e-02 4.954e-01
## [24,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## [25,] 2.368e-01 8.598e-01 6.065e-01 4.482e-27 2.804e-02 3.584e-01
## [26,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## [,13] [,14] [,15] [,16] [,17] [,18]
## [1,] 2.924e+00 1.737e+01 3.945e+00 2.971e+00 3.124e-01 1.595e+01
## [2,] 4.454e-02 5.668e-02 1.223e+00 9.869e-02 7.030e-03 2.071e-01
## [3,] 2.227e-01 8.534e-01 2.930e+00 3.471e-01 2.893e-02 1.134e+00
## [4,] 2.981e-01 9.416e-01 5.042e+00 5.160e-01 4.090e-02 1.483e+00
## [5,] 2.766e+00 1.580e+01 7.311e+00 2.975e+00 3.025e-01 1.498e+01
## [6,] 3.368e-01 1.756e+00 1.830e+00 4.054e-01 3.867e-02 1.795e+00
## [7,] 1.986e-01 8.122e-01 2.324e+00 2.962e-01 2.523e-02 1.020e+00
## [8,] 5.181e-02 -1.378e+00 9.500e+00 4.857e-01 2.395e-02 -7.938e-03
## [9,] 1.378e+00 8.050e+00 2.651e+00 1.437e+00 1.488e-01 7.494e+00
## [10,] -5.424e-24 -3.524e-23 9.507e-24 -4.740e-24 -5.466e-25 -3.011e-23
## [11,] 2.669e-02 1.246e-01 2.259e-01 3.584e-02 3.222e-03 1.397e-01
## [12,] 3.622e-01 1.758e+00 2.694e+00 4.693e-01 4.300e-02 1.907e+00
## [13,] 2.368e-02 -1.461e-01 1.636e+00 9.775e-02 5.663e-03 7.976e-02
## [14,] 1.040e+00 5.689e+00 4.167e+00 1.184e+00 1.165e-01 5.590e+00
## [15,] 2.190e+00 1.305e+01 2.762e+00 2.217e+00 2.336e-01 1.195e+01
## [16,] 7.339e-02 2.184e-01 1.317e+00 1.305e-01 1.022e-02 3.628e-01
## [17,] 5.954e-03 3.507e-02 9.786e-03 6.132e-03 6.395e-04 3.242e-02
## [18,] 5.746e-01 2.273e+00 7.161e+00 8.772e-01 7.385e-02 2.937e+00
## [19,] 3.816e-01 1.647e+00 3.986e+00 5.472e-01 4.754e-02 1.974e+00
## [20,] 4.666e-01 8.286e-01 1.150e+01 9.734e-01 7.107e-02 2.210e+00
## [21,] 4.856e-01 2.891e+00 6.236e-01 4.920e-01 5.181e-02 2.649e+00
## [22,] -1.724e-01 -1.697e+00 3.530e+00 -2.441e-03 -1.107e-02 -1.056e+00
## [23,] -4.820e-02 -9.117e-01 3.433e+00 1.116e-01 1.682e-03 -3.706e-01
## [24,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## [25,] 1.445e-01 7.712e-01 6.859e-01 1.694e-01 1.640e-02 7.733e-01
## [26,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## [,19] [,20] [,21] [,22] [,23] [,24]
## [1,] 7.706e+00 2.226e+01 3.004e+00 6.013e+00 3.105e+00 0
## [2,] 2.125e-01 8.239e-01 1.156e-01 3.942e-03 3.767e-01 0
## [3,] 8.021e-01 2.792e+00 3.868e-01 2.599e-01 9.812e-01 0
## [4,] 1.165e+00 4.203e+00 5.848e-01 2.633e-01 1.630e+00 0
## [5,] 7.582e+00 2.255e+01 3.057e+00 5.418e+00 3.951e+00 0
## [6,] 1.000e+00 3.137e+00 4.286e-01 5.890e-01 7.471e-01 0
## [7,] 6.914e-01 2.368e+00 3.274e-01 2.534e-01 7.931e-01 0
## [8,] 9.075e-01 4.325e+00 6.195e-01 -6.040e-01 2.725e+00 0
## [9,] 3.698e+00 1.082e+01 1.464e+00 2.775e+00 1.688e+00 0
## [10,] -1.292e-23 -3.428e-23 -4.563e-24 -1.242e-23 -9.961e-25 0
## [11,] 8.586e-02 2.823e-01 3.882e-02 4.057e-02 8.212e-02 0
## [12,] 1.135e+00 3.676e+00 5.045e-01 5.786e-01 1.009e+00 0
## [13,] 1.936e-01 8.491e-01 1.207e-01 -7.218e-02 4.795e-01 0
## [14,] 2.968e+00 9.071e+00 1.235e+00 1.931e+00 1.887e+00 0
## [15,] 5.757e+00 1.659e+01 2.239e+00 4.519e+00 2.271e+00 0
## [16,] 2.930e-01 1.066e+00 1.485e-01 5.917e-02 4.228e-01 0
## [17,] 1.584e-02 4.605e-02 6.222e-03 1.211e-02 6.819e-03 0
## [18,] 2.036e+00 7.035e+00 9.738e-01 7.005e-01 2.418e+00 0
## [19,] 1.289e+00 4.351e+00 6.004e-01 5.231e-01 1.388e+00 0
## [20,] 2.119e+00 8.083e+00 1.132e+00 1.403e-01 3.574e+00 0
## [21,] 1.277e+00 3.683e+00 4.969e-01 1.001e+00 5.067e-01 0
## [22,] -1.468e-01 2.559e-01 4.878e-02 -6.380e-01 8.822e-01 0
## [23,] 1.589e-01 1.091e+00 1.605e-01 -3.626e-01 9.393e-01 0
## [24,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0
## [25,] 4.211e-01 1.305e+00 1.780e-01 2.603e-01 2.925e-01 0
## [26,] 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0
## [,25] [,26]
## [1,] 1.508e+00 0
## [2,] 8.485e-02 0
## [3,] 2.548e-01 0
## [4,] 4.005e-01 0
## [5,] 1.617e+00 0
## [6,] 2.469e-01 0
## [7,] 2.118e-01 0
## [8,] 5.283e-01 0
## [9,] 7.531e-01 0
## [10,] -1.903e-24 0
## [11,] 2.387e-02 0
## [12,] 3.041e-01 0
## [13,] 9.755e-02 0
## [14,] 6.836e-01 0
## [15,] 1.119e+00 0
## [16,] 1.026e-01 0
## [17,] 3.164e-03 0
## [18,] 6.360e-01 0
## [19,] 3.814e-01 0
## [20,] 8.189e-01 0
## [21,] 2.488e-01 0
## [22,] 1.112e-01 0
## [23,] 1.612e-01 0
## [24,] 0.000e+00 0
## [25,] 1.007e-01 0
## [26,] 0.000e+00 0Part D
x = (svd(G2)$u)[,2]
y = (svd(G2)$v)[,2]
plot(x,y)
findLetter = function(index){
print(x[index])
print(y[index])
}
findLetter(1) #find A ## [1] 0.509
## [1] -0.1582findLetter(5) #find E ## [1] 0.3304
## [1] -0.5316findLetter(9) #find I## [1] 0.2066
## [1] -0.1791findLetter(15) #find O## [1] 0.3894
## [1] -0.1074findLetter(21) #find U## [1] 0.08585
## [1] 0.05152The vowels are clustered at high-x and low-y values on this plot.
Part E
svdS = svd(S)
sigmaS = spMatrix(26,26,i=seq(1,26),j=seq(1,26), x = rep(0,26))
sigmaS[1,1] = max(svdS$d)
sigmaS[2,2] = sort(svdS$d, partial = 25)[25]
image(sigmaS)
Srank = (svdS$u)%*%sigmaS%*%t(svdS$v)## Note: method with signature 'matrix#Matrix' chosen for function '%*%',
## target signature 'matrix#dgTMatrix'.
## "ANY#TsparseMatrix" would also be validx2 = (svd(Srank)$u)[,2]
y2 = (svd(Srank)$v)[,2]
plot(x2,y2)
This time, the vowels are clustered in the opposite corner (because we can flip the signs on the U and V matrices in the SVD and the relationship between them still holds!)
Problem 4
Here is the vector \(y\) of noisy measurements:
source("http://www.macalester.edu/~dshuman1/data/365/2015.365.Final.Data.r")
print(y)## [1] 40.675 42.337 40.090 41.465 39.768 38.405 39.475 40.928 39.321
## [10] 44.909 39.994 38.778 40.580 38.691 40.404 40.046 37.030 39.838
## [19] 36.835 39.013 39.236 36.457 40.285 38.393 37.054 39.587 39.677
## [28] 38.996 36.966 40.690 36.287 37.034 36.203 35.557 34.159 35.192
## [37] 35.186 36.989 35.092 34.190 33.149 34.989 35.680 33.802 30.487
## [46] 31.636 34.200 31.513 29.230 30.258 32.130 29.711 31.787 31.735
## [55] 27.855 30.170 28.277 29.061 27.366 24.303 24.722 24.631 28.376
## [64] 23.154 24.963 25.875 24.925 23.888 22.781 24.526 23.734 22.241
## [73] 24.427 24.313 19.962 21.316 20.401 19.889 23.061 20.160 18.816
## [82] 18.320 18.445 19.849 17.233 18.246 18.296 16.605 18.375 15.860
## [91] 16.512 14.403 13.254 16.047 15.322 14.619 16.724 16.058 16.186
## [100] 13.637 16.078 14.756 15.450 14.885 13.884 13.181 13.422 13.826
## [109] 13.028 15.453 15.017 12.737 14.074 12.619 11.128 12.291 14.716
## [118] 11.925 12.357 15.901 12.938 12.612 12.940 16.270 12.545 14.012
## [127] 14.545 14.237 14.167 13.139 15.029 14.967 13.785 15.322 16.370
## [136] 15.238 17.360 14.867 14.602 17.665 14.829 15.938 14.806 17.654
## [145] 16.339 14.959 17.467 18.104 18.813 18.201 19.026 22.704 18.822
## [154] 19.753 16.750 20.903 19.045 19.919 23.681 20.723 22.278 24.491
## [163] 21.867 22.389 23.542 22.554 24.408 23.022 24.234 25.882 26.206
## [172] 27.225 23.897 25.401 24.077 26.123 27.653 29.492 25.680 26.963
## [181] 26.516 31.577 27.766 30.263 31.276 29.027 29.720 33.156 29.066
## [190] 31.769 30.161 31.331 30.391 29.570 30.034 31.996 32.648 31.797
## [199] 34.534 33.482 34.396 32.460 31.661 29.980 34.700 33.107 33.584
## [208] 33.641 36.094 36.126 35.367 38.487 34.374 33.614 33.777 35.021
## [217] 34.835 32.298 35.376 37.191 36.305 36.059 34.153 36.209 37.075
## [226] 34.618 34.870 31.757 36.516 33.460 35.565 34.226 32.714 35.839
## [235] 34.247 33.302 30.380 34.941 36.105 33.462 31.645 33.252 35.255
## [244] 32.103 32.690 33.264 32.511 31.503 31.004 32.094 34.109 30.180
## [253] 32.628 33.407 30.454 32.068 28.413 31.369 29.221 29.563 27.947
## [262] 27.766 26.171 25.712 30.644 26.738 25.517 27.924 25.514 23.277
## [271] 25.256 24.939 23.658 27.987 26.172 25.237 24.848 23.342 25.631
## [280] 24.446 22.473 23.962 24.271 25.372 20.571 23.529 23.515 22.258
## [289] 21.036 19.239 20.404 20.087 19.367 20.761 18.753 16.879 21.775
## [298] 19.498 15.062 19.803 18.160 17.217 20.743 17.780 16.421 17.015
## [307] 15.927 16.084 16.610 18.095 15.467 13.843 16.406 16.809 17.053
## [316] 15.793 15.859 13.993 13.762 17.794 15.622 16.240 16.118 14.435
## [325] 17.076 16.730 18.267 18.797 15.705 16.816 16.481 15.992 12.717
## [334] 14.039 15.631 17.382 14.938 14.763 16.116 15.781 13.664 12.548
## [343] 18.473 15.996 15.896 15.799 14.626 13.983 17.747 15.011 14.945
## [352] 16.617 17.807 16.588 17.640 17.376 17.698 19.410 17.339 16.350
## [361] 15.154 19.290 19.198 18.597 19.910 17.633 19.084 17.425 19.252
## [370] 19.885 17.601 17.091 19.051 19.565 21.904 20.509 20.558 22.285
## [379] 21.728 19.481 19.673 20.627 22.597 20.062 23.449 21.977 22.343
## [388] 24.179 25.731 22.936 20.212 22.731 22.693 26.579 22.988 25.652
## [397] 27.606 20.991 26.213 25.902 27.395 22.992 24.590 27.818 25.867
## [406] 23.734 25.634 28.019 23.327 25.587 24.473 26.836 25.269 25.077
## [415] 25.032 23.952 24.650 24.134 27.726 26.674 24.870 25.260 26.069
## [424] 24.431 28.066 26.307 26.416 25.663 27.638 26.092 27.541 27.024
## [433] 27.240 28.748 24.535 25.848 25.194 27.307 24.807 26.344 25.820
## [442] 25.091 24.780 25.347 29.412 28.312 23.111 25.435 24.258 26.435
## [451] 27.267 26.525 23.786 23.838 27.194 26.868 26.512 21.690 27.024
## [460] 23.805 26.158 29.694 22.350 26.097 26.599 23.922 21.778 21.204
## [469] 25.023 24.397 21.469 23.883 24.543 24.329 26.813 25.746 26.616
## [478] 21.872 27.473 23.988 24.264 26.319 23.610 24.183 25.172 23.972
## [487] 25.051 24.955 25.140 25.874 24.294 22.671 21.663 22.555 23.707
## [496] 21.626 24.757 24.582 24.497 25.275 25.406 25.898 24.162 23.790
## [505] 25.000 26.081 23.214 24.889 20.947 26.928 24.741 24.406 26.563
## [514] 23.976 26.801 27.827 26.739 23.466 26.999 26.967 26.912 29.338
## [523] 27.424 29.841 29.248 27.383 28.176 31.207 27.224 26.213 32.298
## [532] 29.561 31.069 31.379 33.774 29.533 29.046 29.520 29.145 27.835
## [541] 28.934 29.119 35.894 33.331 34.467 31.672 31.829 33.106 32.447
## [550] 35.305 36.364 33.666 35.365 33.570 32.643 36.973 33.348 35.206
## [559] 37.078 36.907 38.082 36.648 36.655 38.906 36.938 38.531 36.993
## [568] 38.858 39.301 37.122 41.042 37.116 38.359 35.986 39.553 40.014
## [577] 40.420 39.988 40.068 39.569 38.816 39.699 40.426 40.350 38.781
## [586] 36.294 40.053 39.558 39.659 38.741 37.773 41.730 42.385 40.908
## [595] 41.845 39.408 36.286 40.280 39.199 39.767 42.769 39.576 39.956
## [604] 40.344 40.155 38.558 38.784 38.767 40.700 39.030 36.134 40.227
## [613] 36.690 39.095 36.992 36.983 35.760 37.441 34.621 33.275 35.702
## [622] 36.113 34.529 35.250 35.648 33.397 31.780 33.437 36.303 32.517
## [631] 31.410 35.555 32.953 30.683 33.367 28.043 29.560 27.573 28.317
## [640] 30.171 29.731 27.765 27.248 27.681 27.032 27.974 27.545 24.834
## [649] 28.227 25.040 24.882 21.687 23.818 22.219 23.094 24.879 21.887
## [658] 21.273 22.308 19.325 17.173 17.781 14.784 20.643 16.956 17.483
## [667] 14.937 17.945 15.666 17.767 17.097 16.129 14.274 13.558 15.181
## [676] 12.390 10.458 11.725 12.540 12.522 11.928 12.639 11.975 9.140
## [685] 10.525 8.809 10.126 8.163 10.565 4.931 7.038 5.936 6.590
## [694] 10.438 5.251 7.126 6.062 3.784 5.273 4.937 4.299 4.228
## [703] 5.979 2.974 2.694 5.066 2.238 2.716 4.289 8.114 4.436
## [712] 4.757 4.967 4.913 5.654 2.471 3.317 5.451 3.826 5.352
## [721] 4.131 2.714 5.176 5.179 1.741 3.516 5.873 5.866 4.878
## [730] 9.601 3.696 4.910 7.870 4.456 4.881 7.186 7.275 7.613
## [739] 8.302 7.610 6.911 8.582 8.331 10.579 12.349 9.861 8.983
## [748] 9.300 12.835 11.288 12.569 10.421 11.664 14.514 10.643 15.160
## [757] 14.516 16.625 16.956 15.839 16.884 16.564 15.343 18.206 17.839
## [766] 18.725 19.778 17.161 16.199 19.334 18.263 19.027 20.071 23.701
## [775] 23.201 23.904 23.774 22.710 20.311 24.537 24.703 22.446 25.744
## [784] 27.714 25.356 24.774 25.653 26.425 24.327 27.094 27.593 30.721
## [793] 29.501 28.186 29.201 29.947 30.140 34.234 30.156 30.744 30.222
## [802] 29.644 34.956 32.782 31.900 37.038 32.592 35.406 37.083 33.762
## [811] 36.198 37.653 33.444 34.026 38.667 37.108 39.467 37.012 36.340
## [820] 33.794 34.833 35.042 36.331 37.250 37.896 36.709 37.989 37.474
## [829] 37.272 35.512 39.876 37.834 37.380 35.108 37.737 35.903 38.786
## [838] 35.687 36.439 37.977 36.526 38.071 35.671 35.273 41.028 38.137
## [847] 35.850 36.220 39.384 33.504 35.937 36.514 36.289 35.518 36.949
## [856] 35.133 37.131 32.270 35.542 34.952 32.945 37.526 37.027 33.057
## [865] 34.959 34.515 36.883 33.751 35.185 34.773 34.969 34.367 33.994
## [874] 32.812 32.247 31.409 31.614 32.666 31.182 30.343 29.996 33.012
## [883] 31.097 29.279 28.539 27.625 30.311 28.293 26.161 29.886 28.000
## [892] 26.604 27.031 26.230 27.781 26.376 28.834 27.135 26.886 29.193
## [901] 26.951 25.763 25.155 27.933 26.328 25.629 26.024 25.597 26.672
## [910] 26.613 24.270 24.777 26.645 26.864 24.460 26.291 25.084 23.482
## [919] 22.767 21.826 25.882 22.838 24.970 23.899 23.713 22.558 26.519
## [928] 24.181 23.456 21.947 25.871 26.350 23.106 22.352 21.198 23.125
## [937] 23.666 23.796 25.600 25.770 23.464 24.471 24.345 25.474 23.682
## [946] 25.298 24.903 23.505 24.387 21.400 26.461 21.937 26.749 25.020
## [955] 27.041 26.375 24.323 23.851 24.324 26.972 25.645 25.797 23.342
## [964] 27.316 26.294 26.163 25.752 26.168 28.575 23.923 28.081 23.725
## [973] 27.745 26.749 26.961 25.927 26.207 26.715 26.793 25.329 26.295
## [982] 27.408 29.356 25.191 26.799 25.837 29.618 28.432 27.219 28.899
## [991] 31.108 29.081 25.913 30.472 28.108 31.273 28.040 28.784 29.823
## [1000] 29.628 27.962plot(0:1000,y,type="l",lwd=1,col="blue",xlab="k",ylab="y",ylim=c(0,50))
grid()
Part A
makeA = function(lambda){
d = c(rep(1,1001))
l = c(rep(1,1000))
u = c(rep(0,1000))
Atop = as.matrix(TriDiag(d,l,u))
d2 = c(rep(0, 1001))
l2 = c(rep(1, 1000))
u2 = c(rep(lambda,1000))
Abottom = as.matrix(TriDiag(d2,l2,u2))
A = rbind(Atop,Abottom)
A = A[-2002,]
A
}
A100 = as.matrix(makeA(100))
A10 = as.matrix(makeA(10))
A = as.matrix(makeA(1))
A[1:5,1:5]## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 0
## [2,] 1 1 0 0 0
## [3,] 0 1 1 0 0
## [4,] 0 0 1 1 0
## [5,] 0 0 0 1 1dim(A)## [1] 2001 1001A is a 2001x1001 matrix that accounts for both entries in front of and behind \(y_k\). In that, it puts (lambda) weight on entries behind. (First 5 rows shown)
makeB = function(lambda){
y = as.matrix(y)
y2 = (lambda)*y
b = rbind(y,y2)
b = as.matrix(b[-2002,])
b
}
b100 = as.matrix(makeB(100))
b10 = as.matrix(makeB(10))
b = as.matrix(makeB(10))
b[1:5]## [1] 40.68 42.34 40.09 41.47 39.77dim(b)## [1] 2001 1On the same coin, b is a 1001x1 matrix that accounts for the estimted value of each equation in the A-matrix. Again, the lambda is included to counteract the inclusion of lambda in A. (First 5 rows shown).
Part B — 3 methods
This is essentially a least squares problem because we are trying to fit a line of best fit between each of the two points in the plot. Here are three ways to solve this problem:
Method 1: QR Decomposition: This is the decomposition of the matrix into a product A = QR of an orthogonal matrix ! and an upper triangular matrix R. A must have n linearly indpendent columns to do so.
Method 2: Singular Value Decomposition (SVD): This is useful for solving overdetermined systems of equations. The SVD can be used for computing the pseudoinverse of a matrix A with the singluar value decomposition A = U(Sigma)(V^T)
Method 3: Normal Equations: This is our old, trusty method for solving least squares systems of equations. There are 3 steps to using this method: (1) Choosing a model, (2) Forcing the model to fit the data, and (3) Finally solving the least squares equations.
Part C
solve100 = qr.solve(A100,b100)
solve10 = qr.solve(A10,b10)
solve1 = qr.solve(A,b)
plot(solve100, col = 'yellow')
points(solve10, col = 'blue')
points(solve1, col = 'red')
title(main='Smoothing the Signal')
Problem 5
You can test your function on this matrix \(A\):
A = rbind(c(1,2,3,4), c(2,3,6,5), c(3,6,1,4), c(4,5,4,0))
myEigen = function(A){
eigenval = c() #Empty matrix for eigenvalues
eigenvec = matrix(nrow=nrow(A), ncol = nrow(A)) #Empty matrix for eigenvectors
x = cbind(c(1,0,0,1)) #Arbitrary initial estimate for x
i=1
while(i<(nrow(A)+1)){
eigenval[i] = PI(A,x)$val #Gets eigenvalues
eigenvec[,i] = PI(A,x)$vec #Gets eigenvectors...divide by lamda?
trans = t(PI(A,x)$vec) #middle step for new A
u = (PI(A,x)$vec) #/vnorm(x
right = eigenval[i]*u%*%trans ###
A = A - right #get our new A!
i=i+1 #iterate for new A.
}
return(list(eigenval = eigenval, eigenvec=eigenvec))
}
myEigen(A)## $eigenval
## [1] 13.6036 -4.8892 -3.7465 0.0321
##
## $eigenvec
## [,1] [,2] [,3] [,4]
## [1,] 0.3723 -0.4124 -0.1419 0.8192
## [2,] 0.5945 -0.5643 0.2635 -0.5086
## [3,] 0.5247 0.3688 -0.7454 -0.1819
## [4,] 0.4823 0.6128 0.5957 0.1925eigen(A)## $values
## [1] 13.6036 0.0321 -3.7465 -4.8892
##
## $vectors
## [,1] [,2] [,3] [,4]
## [1,] -0.3723 0.8192 0.1419 0.4124
## [2,] -0.5945 -0.5086 -0.2635 0.5643
## [3,] -0.5247 -0.1819 0.7454 -0.3688
## [4,] -0.4823 0.1925 -0.5957 -0.6128Eigen and myEigen return the same result. Therefore, myEigen works!