Complex-valued linear models

Someone has probably already written code to do this. But I couldn’t find it in CRAN, so here goes.

Oh no lm() won’t take complex numbers! (or rather, it’ll take them, but it’ll discard the imaginary part.)

Easy enough fix. Split into real.

y = a*x
y1 + y2*i = (a1 + a2*i) * (x1 + x2*i)
y1 = a1*x1 - a2*x2
y2 = a1*x2 + a2*x1

So, if we’re looking to discover the coefficients a1 and a2, we can split the vectors y and x like

Re(y1)  Re(x1)  -Im(x1)
Im(y1)  Im(x1)   Re(x1)
...     ...      ...
break1 = function(X) {
        do.call(c, lapply(X, function(x) { c(Re(x), Im(x)) }))
}

break2 = function(X) {
        do.call(c, lapply(X, function(x) { c(-Im(x), Re(x)) }))
}

So if we have a function to do a complex fit, the first thing is to make new variables based on the inputs.

fit.complex = function(Y, X.List) {

        # Split into real variables
        YF = break1(Y)
        XF.List = do.call(c, lapply(X.List,
                function(x) { list(break1(x), break2(x)) } ))

        # ...
}

Then put those into a data.frame and make an appropriate formula.

        # Make the data.fram
        Data = data.frame(Y = YF)
        X.Names = paste('X', 1:length(XF.List), sep='')

        for (N in seq_along(XF.List)) {
                Data[[ X.Names[[N]] ]] = XF.List[[N]]
        }

        Formula = paste("Y ~ ", paste(X.Names, collapse='+'), "-1")

It’s important to put the “-1” in the formula so that lm() doesn’t include a constant term. (We might want a constant term, but it would have to look more like c(1, 0, 1, 0, …) because it is complex).

Then do the fit and extract the coefficients.

        Model = lm(as.formula(Formula), data=Data)

        # Make them complex again
        Coeffs = sapply(seq_along(X.List),
                function(N) {
                        ( Model$coefficients[[ X.Names[[2*N-1]] ]]
                        + Model$coefficients[[ X.Names[[2*N]] ]]*1i )
                })
        names(Coeffs) = names(X.List)

        Model$coefficients.complex = Coeffs

The whole function looks like

fit.complex = function(Y, X.List) {

        # Split into real variables
        YF = break1(Y)
        XF.List = do.call(c, lapply(X.List,
                function(x) { list(break1(x), break2(x)) } ))

        # Make the data.fram
        Data = data.frame(Y = YF)
        X.Names = paste('X', 1:length(XF.List), sep='')

        for (N in seq_along(XF.List)) {
                Data[[ X.Names[[N]] ]] = XF.List[[N]]
        }

        # Formula + Model
        Formula = paste("Y ~ ", paste(X.Names, collapse='+'), "-1")
        Model = lm(as.formula(Formula), data=Data)

        # Make them complex again
        Coeffs = sapply(seq_along(X.List),
                function(N) {
                        ( Model$coefficients[[ X.Names[[2*N-1]] ]]
                        + Model$coefficients[[ X.Names[[2*N]] ]]*1i )
                })
        names(Coeffs) = names(X.List)

        Model$coefficients.complex = Coeffs

        Model
}

Now test it

Beta0 = 1 + 3i
Beta1 = 3 - 2i

X = runif(15, 0, 10)
Y = (Beta0 + Beta1*X +
        rnorm(length(X), 0, 0.7) * exp(1i*runif(length(X), 0, 2*pi))
)

Model = fit.complex(Y, list(
         const = 0*X+1,
        linear = X
))

Beta0.Est = Model$coefficients.complex[[1]]
Beta1.Est = Model$coefficients.complex[[2]]
> Beta0.Est
[1] 1.090385+3.017922i
> Beta1.Est
[1] 2.912617-2.030427i

Excellent.

Advertisements
This entry was posted in R. Bookmark the permalink.

2 Responses to Complex-valued linear models

  1. Reygun says:

    great”!

  2. hello!,I love your writing very so much! percentage we communicate extra about your post on AOL? I need an expert in this area to solve my problem. Maybe that’s you! Having a look forward to peer you.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s