[R] Regression of complex-valued functions
Andrea Graziani
a.graziani at mail.com
Thu Feb 13 00:25:53 CET 2014
Hi Frede,
Thank you for your accurate answer!
If I understand well, your way to use nls() solves the problem using too many physical parameters.
I solved the problem following the other way that you and Rolf Turner suggested (i.e. splitting the complex-valued problem into two real-valued problems and using nlm() ).
If you are interested in more details, please see my answer to Rolf’s mail.
> 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?
You noticed a key point of the mechanical problem!
However I’m looking for a “single curve”, so-called “master curve”, to fit all experimental data!
>
> 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).
;)
Thanks again!
Andrea
>
> ## 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.
> <plot-raw-data.png><fitted-Re.png><fitted-abs.png>
More information about the R-help
mailing list