Regression of complex-valued functions

Hi Andrea

Although there are theory on linear models for complex valued normal distributions (see e.g. http://www.amazon.com/Linear-Graphical-Models-Multivariate-Distribution/dp/0387945210) to my knowledge there is no function for doing complex-valued regression in R.

However as my reference probably explain (long time since I saw the textbook) you can consider the Im and Re parts separately. But first some questions. I don't know anything about the physics behind your formula. It seems that you want to do non-linear regression for five curves (one for each value of components of the vector logaT). Why are loga40, loga30, loga20, 0 (zero value), and loga0 not independent variables that you control? If so you don't want to estimate those.

If the logaT parameters are not under your control do you still really want to estimate the same set of values of the parameters of the model for each set of your data? And not separate estimates for each curve?

However (putting my head under my shoulder) you can do something like this (Having just noticed Rolf Turner's answer you should absolutely explore that solutions. That would probably the best way to go).

###### R script:

## added some NAs do make separat plot of curves
## then you need to set na.action argument of nls()

## your data
f_data <- rep(c(0.1,0.25,1,4,12, NA),5)

E <- c(289.7411+ 225.0708i , 386.4417+ 303.5021i ,  671.5132+ 521.1253i , 1210.8638+ 846.6253i ,     1860.9623+1155.4139i,NA,
       862.8984+ 636.2637i , 1159.0436+ 814.5609i , 1919.0369+1186.5679i , 3060.7207+1573.6088i  ,   4318.1781+1868.4761i,NA,
       2760.7782+1418.5450i , 3306.3013+1612.2712i , 4746.6958+1923.8662i , 6468.5769+2148.9502i ,   8072.2642+2198.5344i,NA,
       6757.7680+2061.3110i , 7591.9787+2123.9168i , 9522.9471+2261.8489i , 11255.0952+2166.6411i , 12601.3970+2120.7178i,NA,
       11913.6543+2016.0828i , 12906.8294+2030.0610i ,14343.7693+1893.4877i, 15942.7703+1788.0910i, 16943.2261+1665.9847i,NA)

## function
E_2S2P1D <- function(f,logaT,Eg,Ee,k,h,delta,logbeta,logtau) {
    Ee+(Eg-Ee)*(
        1+delta*(2i*pi*10^(logtau)*f*10^logaT)^-k +
        (2i*pi*10^(logtau)*f*10^logaT)^-h +
        (2i*pi*10^(logtau)*f*10^logaT*10^(logbeta))^-1
        )^-1 }

mycol <- rep(c("red", "green", "blue", "magenta", "black"), each = 6)

png("plot-raw-data.png", height = 3*480)
par(mfrow = c(3,1))
plot(f_data, Re(E), type = "p", col =  mycol)# log="xy")
plot(f_data,Im(E), type = "p", col =  mycol)# log="xy")
plot(f_data,abs(E), type = "p", col = mycol)# log="xy")
dev.off()

myE <- data.frame(f = f_data, E = E)

## this seems to run
fit.re <- nls(Re(E) ~ Re(E_2S2P1D(f,c(rep(loga40,5),rep(loga30,5),rep(loga20,5),rep(0,5),rep(loga0,5)),
                                      Eg,Ee,k,h,delta,logbeta,logtau)),
              data = myE, na.action = "na.exclude",  
              start=list(loga40=-3.76,loga30=-2.63,loga20=-1.39,loga0=1.68,
                  Eg=27000,Ee=200,k=0.16,h=0.47,delta=1.97,logbeta=2.2,logtau=-0.02) )

fit.re

## plot fit overlayed observations
png("fitted-Re.png")
plot(f_data, Re(E), type = "p", col = mycol)# log="xy")
lines(f_data, fitted(fit.re), col = "grey80") #, log = "xy"
dev.off()

## this gives an error
## probably due to the 5 curves do not share parameters. See 3rd panel in attached plot-raw-data.png
fit.im <- nls(Im(E) ~ Im(E_2S2P1D(f,c(rep(loga40,5),rep(loga30,5),rep(loga20,5),rep(0,5),rep(loga0,5)),
                                  Eg,Ee,k,h,delta,logbeta,logtau)),
              data = myE, na.action = "na.exclude",
              start=list(loga40=-3.76,loga30=-2.63,loga20=-1.39,loga0=1.68,
                   Eg=27000,Ee=200,k=0.16,h=0.47,delta=1.97,logbeta=2.2,logtau=-0.02) )

## Now try abs()
## which works and gives allmost same estimates as for fit.re
fit.abs <- nls(abs(E) ~ abs(E_2S2P1D(f,c(rep(loga40,5),rep(loga30,5),rep(loga20,5),rep(0,5),rep(loga0,5)),
                                      Eg,Ee,k,h,delta,logbeta,logtau)),
              data = myE, na.action = "na.exclude",
              start=list(loga40=-3.76,loga30=-2.63,loga20=-1.39,loga0=1.68,
                  Eg=27000,Ee=200,k=0.16,h=0.47,delta=1.97,logbeta=2.2,logtau=-0.02))

fit.abs

png("fitted-abs.png")
plot(f_data, abs(E), type = "p", col = mycol)# log="xy")
lines(f_data, predict(fit.abs), col = "grey80") #, log = "xy"
dev.off()

Yours sincerely / Med venlig hilsen

Frede Aakmann T?gersen
Specialist, M.Sc., Ph.D.
Plant Performance & Modeling

Technology & Service Solutions
T +45 9730 5135
M +45 2547 6050
frtog at vestas.com
http://www.vestas.com

Company reg. name: Vestas Wind Systems A/S
This e-mail is subject to our e-mail disclaimer statement.
Please refer to www.vestas.com/legal/notice
If you have received this e-mail in error please contact the sender.

Regression of complex-valued functions

Thread (8 messages)