# Complex-valued linear models

September 16, 2010
By

(This article was first published on Struggling Through Problems » R, and kindly contributed to R-bloggers)

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[]
Beta1.Est = Model\$coefficients.complex[]```
```> Beta0.Est
 1.090385+3.017922i
> Beta1.Est
 2.912617-2.030427i```

Excellent.

To leave a comment for the author, please follow the link and comment on their blog: Struggling Through Problems » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...

If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Tags:

Comments are closed.

# Never miss an update! Subscribe to R-bloggers to receive e-mails with the latest R posts.(You will not see this message again.)

Click here to close (This popup will not appear again)