4

I've gotten a MLE estimator that is written in GAUSS that I am trying to re-code into R. I do not use, nor have I ever used GAUSS itself (and do not have access to it). In the code, there is a line that is a bit confusing to me.

In the validated GAUSS code there is a line under the "inputs" (part of the comments) that states:

invsig: scalar or m-by-m matrix with inverse of sigma

I am working on getting the code to work piece by piece, but my first question is something relatively simple.

Here is the GAUSS snippet that confuses me:

...
local m, k, tobs, invsig
m = rows(y);k = rows(x); tobs = rows(dat)
invsig= eye(m)*invsig
...

I understand that this is the identity matrix multiplied by the "input" invsig, but in simulated examples that works (from log files attached to the code), the program can be started using the scalar value of invsig. IE: setting an initial value as invsig = 1

In R, this is not working. Here is the simple "test" code to try to get this:

y.mat <- rep(rexp(3)) 
x.mat <- matrix(rexp(36), 12, 3) 
myfct <- function(x,invsig){
   m <- nrow(x)
   invsig <- diag(m)%*%invsig
   return(invsig)
}
t1 <- myfct(x.mat, 1)     ##Non-conformable error
t2 <- myfct(x.mat, y.mat) ##Works

I understand the non-conformable error that I am getting in R. The question is am I missing something in the conversion between GAUSS and R? In reading the help manuals online, GAUSS does matrix operations by using the individual symbols (*/+-) and to do things element-wise, you add a "." before each operation. So to me, the GAUSS code is saying to do matrix multiplication (%*% in R)and that is what the simple function tries to do.

Any comments or suggestions are greatly appreciated!

Bill the Lizard
  • 398,270
  • 210
  • 566
  • 880
Tony
  • 149
  • 1
  • 9
  • This isn't the problem you ask about, but in GAUSS, it looks like `m` is defined as the number of rows in `y`, but in R you make `m` the number of rows in `x`. – Gregor Thomas May 28 '13 at 22:29
  • Thanks, yeah, the sample function is an adaptation of the GAUSS code, (not word for word replica). They have several different matrices like this in a large program, so I accidentally mixed variables, but the problem remains if I go back and change it back from x to y. – Tony May 28 '13 at 22:50

1 Answers1

3

Your problem is that you can't matrix multiply a 3x3 matrix with a 1x1 matrix. I'd recommend something like

myfct <- function(x,invsig){
    if (is.matrix(invsig)) return(invsig)
    m <- nrow(x)
    return(diag(invsig, nrow = m))
}

Edited I had left in the multiplying by diag(m) in the first case, but of course it's unnecessary if invsig is already a matrix.

Gregor Thomas
  • 136,190
  • 20
  • 167
  • 294
  • Thanks I will try this. I am still confused on the difference between GAUSS and R in their matrix multiplication. I don't understand how in their GAUSS code it can both work as a matrix and a scalar at the same time unless I am omitting something. The fact that it can be both a scalar and a length(m) matrix and still solve both ways is unusual. – Tony May 28 '13 at 22:53
  • I don't know GAUSS either, but maybe they have a "nice" and confusing default where a matrix times (`%*%`) a scalar returns the matrix regular times the scalar. – Gregor Thomas May 28 '13 at 22:58
  • GAUSS's approach does actually make mathematical sense; it would make sense to overload `*` as meaning elementwise multiplication in the case of (scalar*matrix) and standard multiplication in the case of (matrix*matrix), because elementwise multiplication *is* the standard mathematical interpretation of scalar-by-matrix multiplication. MATLAB, GAUSS, AD Model Builder are seem more consistent with math notation to me, even if R's Hadamard product notation is more convenient for programming ... – Ben Bolker May 29 '13 at 00:02
  • @BenBolker agreed. Although, definitely something I can see both sides to, like the default `[` of `drop = TRUE`. – Gregor Thomas May 29 '13 at 06:05
  • @BenBolker: Thanks for the clarification. While it might make sense mathematically, it makes things a bit more complex to go between the two languages. @ shujaa: the method works great to work between the scalar/matrix problem. I appreciate it! – Tony May 29 '13 at 11:17