CODEWITHSUNDEEP
r
user 31466
user 31466

Reputation: 699

Non-conformable arrays in R

y <- matrix(c(7, 9, -5, 0, 2, 6), ncol = 1)
try <- t(y)
tryy <- try %*% y
i <- solve(tryy)
h <- y %*% i %*% try

uniroot(as.vector(solve(((1-x) * diag(6)) + h)), c(-Inf, Inf))

Error in (1 - x) * diag(6) : non-conformable arrays

The purpose of this command uniroot(as.vector(solve(((1-x) * diag(6)) + h)), c(-Inf, Inf)) is to solve the characteristics equation det[(1-λ)I+h] = 0

where, λ=eigenvalues , I=identity matrix , h=hat matrix=y(y'y)^(-1)y'

here λ is unknown ,we have to solve for it.

I am not understanding where is the problem here? I have tried as:

  as.vector(solve(6*diag(6)+h))

This is not non-conformable. But why is not working inside the uniroot function?

Upvotes: 2

Views: 5109

Answers (1)

jlhoward
jlhoward

Reputation: 59355

Your question is a bit confusing, so I have to make a couple of assumptions. If you want the eigenvalues of h, then the characteristic equation is:

det(h - I*λ) = 0

not

det[(1-λ)I+h] = 0

So I used the former.

Given the above, the short answer is: do it this way.

f <- function(lambda) det(h -lambda*diag(6))
F <- Vectorize(f)                         
library(rootSolve)
uniroot.all(F,c(-1000,1000),n=2000)
# [1] 0 1

# or, much more simply
eigen(h)$values
# [1]  1.000000e+00  2.220446e-16  0.000000e+00 -2.731318e-18 -6.876381e-18 -7.365903e-17

So h has 2 eigenvalues, 0 and 1. Note that the built-in function eigen(...) finds 6 roots, but 5 of them are within the machine tolerance of 0.

The question about why your code fails is a bit more involved.

First, your code:

tryy <- try %*% y

is the dot product of y with itself (so, a scalar), returned as a matrix with one element. When you "invert" that using solve(...)

i <- solve(tryy)

you simply take the reciprocal, so i is also a matrix with 1 element. I'm not sure if this is what you had in mind.

Second, uniroot(...) does not work this way. The first argument must be a function; you've passed an expression which depends on x, which in turn is undefined. You could try:

f <- function(x) det(h-x*diag(6))
uniroot(f,c(-Inf,Inf))

but this wouldn't work either because (a) uniroot(...) works on a finite interval, (b) it requires that the function f(...) have different sign at the ends of the interval, and (c) in any event it would return only one root (the smaller one).

So you could use uniroot.all(...) in package rootSolve. uniroot.all(...) also requires a function as it's first argument, but there's a twist: the function must be "vectorized". This means that if you pass a vector of lambda values, f(...) should return a vector of the same length. Fortunately in R there is an easy way to "vectorize" a given function, as in:

F <- Vectorize(f).

Even this has it's limits. uniroot.all(...) also requires a finite interval, so we have to guess what that is, and also it evaluates F on n sub-intervals. So if your interval does not contain all the roots, or if the sub-intervals are not small enough, you will not find all the roots.

Using the built-in eigen(...) function is definitely the best option.

Upvotes: 2

Related Questions