[R] Regression of complex-valued functions

Frede Aakmann Tøgersen frtog at vestas.com
Tue Feb 11 09:50:47 CET 2014


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. 


> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
> On Behalf Of Andrea Graziani
> Sent: 9. februar 2014 23:46
> To: r-help at r-project.org
> Subject: [R] Regression of complex-valued functions
> 
> Hi everyone,
> 
> I previously posted this question but my message was not well written and
> did not contain any code so I will try to do a better job this time.
> 
> The goal is to perform a non-linear regression on complex-valued data.
> I will first give a short description of the data and then describe the complex-
> valued non-linear function.
> 
> ***The Data
> Obtained from mechanical tests performed at 5 loading frequencies 0.1, 0.25,
> 1, 0.4 and 12 Hz, and at 5 temeperatures.
> The independent variable used in the regression is:
> f_data <- rep(c(0.1,0.25,1,4,12),5)
> 
> The measured values of the response variable are:
> E <- c(289.7411+ 225.0708i , 386.4417+ 303.5021i ,  671.5132+ 521.1253i ,
> 1210.8638+ 846.6253i , 1860.9623+1155.4139i ,
> 862.8984+ 636.2637i , 1159.0436+ 814.5609i , 1919.0369+1186.5679i ,
> 3060.7207+1573.6088i  , 4318.1781+1868.4761i  ,
> 2760.7782+1418.5450i , 3306.3013+1612.2712i , 4746.6958+1923.8662i ,
> 6468.5769+2148.9502i , 8072.2642+2198.5344i  ,
> 6757.7680+2061.3110i , 7591.9787+2123.9168i , 9522.9471+2261.8489i ,
> 11255.0952+2166.6411i , 12601.3970+2120.7178i ,
> 11913.6543+2016.0828i , 12906.8294+2030.0610i , 14343.7693+1893.4877i ,
> 15942.7703+1788.0910i , 16943.2261+1665.9847i)
> 
> To visualize the data:
> plot(f_data,Re(E),log="xy")
> plot(f_data,Im(E),log="xy")
> plot(E)
> 
> ***Non-linear regression function:
> Obtained from an analytical model
> 
> 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
> 
> E_2S2P1D is a complex-valued function (note the imaginary unit "i") where:
> "f" is the real-valued independent variable (i.e. the frequency),
> "logaT","Eg","Ee","k","h","delta","logbeta","logtau" are real-valued
> parameters.
> 
> For example:
> 
> > E_2S2P1D(1,0,27000,200,0.16,0.47,1.97,2.2,0.02)
> [1] 9544.759+2204.974i
> 
> > E_2S2P1D(1,-1,27000,200,0.16,0.47,1.97,2.2,0.02)
> [1] 6283.748+2088.473i
> 
> In order to find the parameters of the regression function I use "nls".
> 
> Executing
> 
> fit <- nls(E ~
> E_2S2P1D(f_data,c(rep(loga40,5),rep(loga30,5),rep(loga20,5),rep(0,5),rep(lo
> ga0,5)),
>                             Eg,Ee,k,h,delta,logbeta,logtau),
>           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) )
> 
> results in a lot of warnings (my translation into english):
> "In numericDeriv(form[[3L]], names(ind), env) :
> imaginary parts removed during the conversion"
> 
> I've the same error trying:
> y <- E_2S2P1D(f_data,c(rep(-3.76,5),rep(-2.63,5),rep(-
> 1.39,5),rep(0,5),rep(1.68,5)),27000,200,0.16,0.47,1.97,2.2,0.02)
> plot(E,y)
> 
> If I got it right R is not able to handle regression problems on complex-valued
> functions, correct?
> If yes, is there some workaround?
> For example using MSExcel I write the Real and Imaginary parts separately,
> and than I minimize
> 
> error = sum [ ( Re(E) - Re(E_E_2S2P1D) ) / Re(E) )^2 ] + sum [ ( Im(E) -
> Im(E_E_2S2P1D) ) / Im(E) )^2 ]
> 
> using the MSExcel solver function.
> 
> thank you for your help
> Andrea
> 
> PS
> Thanks to David Winsemius for his suggestions
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
> and provide commented, minimal, self-contained, reproducible code.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: plot-raw-data.png
Type: image/png
Size: 22434 bytes
Desc: plot-raw-data.png
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20140211/1e0e282e/attachment-0006.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fitted-Re.png
Type: image/png
Size: 15021 bytes
Desc: fitted-Re.png
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20140211/1e0e282e/attachment-0007.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fitted-abs.png
Type: image/png
Size: 15622 bytes
Desc: fitted-abs.png
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20140211/1e0e282e/attachment-0008.png>


More information about the R-help mailing list