[R] Regression of complex-valued functions

[David Winsemius]

On Feb 9, 2014, at 2:45 PM, Andrea Graziani wrote:

> 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 

Others have pointed out the theoretical difficulties regarding curve fitting on the complex plane. My concern looking at this is that you have apparently lapse into mathematical notation rather than functional programming notation above:

This 
2i*pi*10^(logtau)*f*10^logaT

Would need to be this:

 2*i*pi*10^(logtau)*f*10^logaT

You would start getting parser errors since '2i' is not a valid token name. I would also worry about negative powers. They can sometimes be numerically unstable. 

And you have not offered the values for temperatures which I suspect are the items in your start list: loga40=-3.76,loga30=-2.63,loga20=-1.39,loga0=1.68. I would think these are fixed rather than something to be optimized. You can leave them in the global environment or build them into the function but I don't think they should be in your start list. (Frede Aakmann Tøgersen noted that issue but did not build a solution that incorporated addressed that concern.)

To examine that relationship you would:

loga40=-3.76;loga30=-2.63;loga20=-1.39;loga0=1.68
logaT <- c(rep(loga40,5),rep(loga30,5),rep(loga20,5),rep(0,5),rep(loga0,5))
 # To visualize the data:

 plot(logaT,Re(E), log="y")
 plot(logaT,Im(E), log="y")

-- 
David.


> 
> 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(loga0,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.

David Winsemius
Alameda, CA, USA