2
$\begingroup$

I have collected this empirical data, and would like to know which function describes this data (3 curves) at best. I have try non linear fit in R, but not good fit. How can I discover the relation between y and x in R?

enter image description here

The data:

"n";"d1";"d2";"d3"

3;0.517638090205041;1;1.4142135623731
4;0.390180644032256;0.765366864730179;1.1111404660392
5;0.312868930080462;0.618033988749895;0.907980999479094
6;0.261052384440103;0.517638090205041;0.76536686473018
7;0.223928952206616;0.445041867912629;0.660558123910334
8;0.196034280659121;0.390180644032257;0.580569354508924
9;0.174311485495316;0.34729635533386;0.517638090205041
10;0.15691819145569;0.312868930080462;0.466890727711811
11;0.142678366398465;0.284629676546571;0.425130579105954
12;0.130806258460286;0.261052384440103;0.390180644032256
13;0.120756994844572;0.241073360510646;0.360510075627812
14;0.112140894474383;0.223928952206615;0.335012446609472
15;0.104671912485888;0.209056926535305;0.31286893008046
16;0.0981353486548355;0.196034280659121;0.293460948910723
17;0.092366917291479;0.184536718926604;0.276312709903764
18;0.0872387747306717;0.174311485495316;0.261052384440102
19;0.0826499484976261;0.165158690944665;0.247385262538695
20;0.0785196315181367;0.156918191455691;0.235074794915676
21;0.0747823885526514;0.149460187172849;0.223928952206616
22;0.0713846676779598;0.142678366398463;0.213790243130225
23;0.0682822203719363;0.136484826729343;0.204528297884069
24;0.065438165643552;0.130806258460286;0.196034280659122
25;0.0628215181562557;0.125581039058627;0.188216626637028
26;0.0604060556017779;0.120756994844572;0.180997751659275
27;0.0581694374862225;0.116289657820951;0.174311485495316
28;0.056092512551737;0.112140894474384;0.168101049858496
29;0.0541587693522699;0.108277817170835;0.162317451054862
30;0.0523538966157466;0.104671912485888;0.15691819145569
31;0.0506654286263763;0.101298337677426;0.151866228450493
32;0.0490824570458241;0.0981353486548363;0.147129127199335
33;0.0475953950922126;0.095163831647483;0.142678366398463
34;0.0461957833060012;0.0923669172914789;0.138488765788384
35;0.0448761285916098;0.0897296607010294;0.134538011224079
36;0.0436297700691224;0.0872387747306735;0.130806258460286
37;0.0424507666672487;0.0848824063922969;0.127275800850932
38;0.0413338024551093;0.0826499484976292;0.123930789257192
39;0.0402741065308789;0.0805318802188283;0.120756994844571
40;0.0392673849212574;0.0785196315181373;0.117741607302378
41;0.0383097624422297;0.0766054673800711;0.114873062467447
42;0.0373977328618061;0.0747823885526547;0.112140894474385
43;0.0365281160138499;0.0730440461153178;0.109535608457286
44;0.0356980207576916;0.0713846676779614;0.107048570547228
45;0.0349048128745667;0.0697989934050007;0.104671912485888
46;0.0341460871499514;0.0682822203719398;0.10239844863967
47;0.0334196430186013;0.0668299540153464;0.100221603574238
48;0.0327234632529732;0.0654381656435497;0.0981353486548354
49;0.032055695260346;0.06410315514331;0.0961341463872389
50;0.0314146346236404;0.062821518156254;0.0942129014192835
51;0.0307987105780134;0.0615901171123406;0.0923669172914777
52;0.0302064731627722;0.0604060556017765;0.0905918581651677
53;0.0296365818273866;0.0592666556451206;0.0888837148726716
54;0.0290877953031653;0.0581694374862237;0.0872387747306751
55;0.0285589625795756;0.0571121015873938;0.0856535946392055
56;0.0280490148471273;0.0560925125517377;0.0841249770565297
57;0.0275569582881908;0.0551086847363239;0.0826499484976255
58;0.0270818676134442;0.0541587693522675;0.0812257402524248
59;0.0266228802554294;0.0532410428755501;0.0798497710610369
60;0.0261791911426893;0.0523538966157475;0.078519631518138
61;0.0257500479877963;0.0514958273099806;0.0772330700085603
62;0.0253347470312821;0.0506654286263776;0.0759879800014808
63;0.0249326291908296;0.0498613834761447;0.0747823885526511
64;0.0245430765714384;0.0490824570458235;0.0736144458827192
65;0.0241655092976361;0.0483274904722663;0.0724824159159856
66;0.0237993826336305;0.047595395092214;0.0713846676779606
67;0.0234441843613305;0.0468851472065193;0.070319667462093
68;0.023099432389684;0.0461957833060036;0.0692859716866716
69;0.0227646725718074;0.0455263957122102;0.0682822203719357
70;0.0224394767092889;0.04487612859161;0.0673071311755866
71;0.0221234407250315;0.0442441743063726;0.0663594939315987
72;0.0218161829883637;0.0436297700691228;0.0654381656435512
73;0.0215173427778236;0.0430321948724444;0.0645420658889694
74;0.0212265788685137;0.0424507666672445;0.0636701725957649
75;0.0209435682324906;0.0418848397667129;0.0628215181562556
76;0.020668004841735;0.0413338024551121;0.0619951858474861
77;0.0203995985643589;0.040797074782811;0.0611903065302504
78;0.0201380741457537;0.0402741065308831;0.0604060556017787
79;0.0198831702670611;0.0397643753301354;0.059641650179499
80;0.0196346386743029;0.0392673849212618;0.0588963464959308
81;0.019392243371958;0.0387826635436474;0.0581694374862219
82;0.0191557598755964;0.038309762442233;0.0574602505521723
83;0.0189249745184607;0.0378482544820363;0.0567681454875653
84;0.018699683807684;0.037397732861809;0.0560925125517379
85;0.0184796938257785;0.036957809918259;0.0554327706788301
86;0.0182648196740942;0.0365281160138496;0.0547883658119943
87;0.0180548849544306;0.036108298501119;0.054158769352271
88;0.0178497212861548;0.0356980207576873;0.0535434767132445
89;0.0176491678558172;0.0352969612862607;0.0529420059729848
90;0.017453070996745;0.0349048128745662;0.0523538966157477
91;0.0172612837964412;0.0345212818107782;0.0517787083566959
92;0.0170736657294776;0.0341460871499508;0.0512160200430905
93;0.0168900823141926;0.0337789600278608;0.0506654286263748
94;0.0167104047912038;0.0334196430186019;0.050126548199817
95;0.0165345098222458;0.0330678895327165;0.0495990090968606
96;0.0163622792078739;0.0327234632529711;0.0490824570458204
97;0.016193599622553;0.0323861376050035;0.0485765523768222
98;0.0160283623660578;0.0320556952603521;0.0480809692772387
99;0.0158664631298188;0.0317319276696139;0.0475953950922138
100;0.0157078017774206;0.0314146346236388;0.0471195296672208
$\endgroup$
1
  • $\begingroup$ try to plot 1/d vs n for each of d1,d2 and d3. You might like to consider that the slopes appear to be all related to $\pi$. $\endgroup$ – Glen_b Sep 24 '17 at 1:11
2
$\begingroup$

The functions that you are seeking are of the form

$y=2 \, \mathrm{ sin} \left(\frac{a \pi}{x} \right)$.

And you got a=1/4, a=2/4 and a=3/4.

(I found this function by realizing that your values resemble much the length of sides of polygons whose vertices are on a unit circle)

The nls function (I assume this is what you tried as 'non linear fit in R') won't work well on your data since it seems to be generated based on calculations and does not include an error term (the nls function uses this error term in approximations of convergence and decisions to terminate).

This issue has been mentioned several times on this website. See for instance: Fiting non-linear model to small data set

Using your data as an example we can show three solutions.

  • use the optim package
  • force the nls function to give it's output even if it has an error
  • add a tiny bit of random noise to the dependent variable

code:

# data set

dv<- c(3,0.517638090205041,1,1.4142135623731,4,0.390180644032256,0.765366864730179,1.1111404660392,5,0.312868930080462,0.618033988749895,0.907980999479094,6,0.261052384440103,0.517638090205041,0.76536686473018,7,0.223928952206616,0.445041867912629,0.660558123910334,8,0.196034280659121,0.390180644032257,0.580569354508924,9,0.174311485495316,0.34729635533386,0.517638090205041,10,0.15691819145569,0.312868930080462,0.466890727711811,11,0.142678366398465,0.284629676546571,0.425130579105954,12,0.130806258460286,0.261052384440103,0.390180644032256,13,0.120756994844572,0.241073360510646,0.360510075627812,14,0.112140894474383,0.223928952206615,0.335012446609472,15,0.104671912485888,0.209056926535305,0.31286893008046,16,0.0981353486548355,0.196034280659121,0.293460948910723,17,0.092366917291479,0.184536718926604,0.276312709903764,18,0.0872387747306717,0.174311485495316,0.261052384440102,19,0.0826499484976261,0.165158690944665,0.247385262538695,20,0.0785196315181367,0.156918191455691,0.235074794915676,21,0.0747823885526514,0.149460187172849,0.223928952206616,22,0.0713846676779598,0.142678366398463,0.213790243130225,23,0.0682822203719363,0.136484826729343,0.204528297884069,24,0.065438165643552,0.130806258460286,0.196034280659122,25,0.0628215181562557,0.125581039058627,0.188216626637028,26,0.0604060556017779,0.120756994844572,0.180997751659275,27,0.0581694374862225,0.116289657820951,0.174311485495316,28,0.056092512551737,0.112140894474384,0.168101049858496,29,0.0541587693522699,0.108277817170835,0.162317451054862,30,0.0523538966157466,0.104671912485888,0.15691819145569,31,0.0506654286263763,0.101298337677426,0.151866228450493,32,0.0490824570458241,0.0981353486548363,0.147129127199335,33,0.0475953950922126,0.095163831647483,0.142678366398463,34,0.0461957833060012,0.0923669172914789,0.138488765788384,35,0.0448761285916098,0.0897296607010294,0.134538011224079,36,0.0436297700691224,0.0872387747306735,0.130806258460286,37,0.0424507666672487,0.0848824063922969,0.127275800850932,38,0.0413338024551093,0.0826499484976292,0.123930789257192,39,0.0402741065308789,0.0805318802188283,0.120756994844571,40,0.0392673849212574,0.0785196315181373,0.117741607302378,41,0.0383097624422297,0.0766054673800711,0.114873062467447,42,0.0373977328618061,0.0747823885526547,0.112140894474385,43,0.0365281160138499,0.0730440461153178,0.109535608457286,44,0.0356980207576916,0.0713846676779614,0.107048570547228,45,0.0349048128745667,0.0697989934050007,0.104671912485888,46,0.0341460871499514,0.0682822203719398,0.10239844863967,47,0.0334196430186013,0.0668299540153464,0.100221603574238,48,0.0327234632529732,0.0654381656435497,0.0981353486548354,49,0.032055695260346,0.06410315514331,0.0961341463872389,
       50,0.0314146346236404,0.062821518156254,0.0942129014192835,51,0.0307987105780134,0.0615901171123406,0.0923669172914777,52,0.0302064731627722,0.0604060556017765,0.0905918581651677,53,0.0296365818273866,0.0592666556451206,0.0888837148726716,54,0.0290877953031653,0.0581694374862237,0.0872387747306751,55,0.0285589625795756,0.0571121015873938,0.0856535946392055,56,0.0280490148471273,0.0560925125517377,0.0841249770565297,57,0.0275569582881908,0.0551086847363239,0.0826499484976255,58,0.0270818676134442,0.0541587693522675,0.0812257402524248,59,0.0266228802554294,0.0532410428755501,0.0798497710610369,60,0.0261791911426893,0.0523538966157475,0.078519631518138,61,0.0257500479877963,0.0514958273099806,0.0772330700085603,62,0.0253347470312821,0.0506654286263776,0.0759879800014808,63,0.0249326291908296,0.0498613834761447,0.0747823885526511,64,0.0245430765714384,0.0490824570458235,0.0736144458827192,65,0.0241655092976361,0.0483274904722663,0.0724824159159856,66,0.0237993826336305,0.047595395092214,0.0713846676779606,67,0.0234441843613305,0.0468851472065193,0.070319667462093,68,0.023099432389684,0.0461957833060036,0.0692859716866716,69,0.0227646725718074,0.0455263957122102,0.0682822203719357,70,0.0224394767092889,0.04487612859161,0.0673071311755866,71,0.0221234407250315,0.0442441743063726,0.0663594939315987,72,0.0218161829883637,0.0436297700691228,0.0654381656435512,73,0.0215173427778236,0.0430321948724444,0.0645420658889694,74,0.0212265788685137,0.0424507666672445,0.0636701725957649,75,0.0209435682324906,0.0418848397667129,0.0628215181562556,76,0.020668004841735,0.0413338024551121,0.0619951858474861,77,0.0203995985643589,0.040797074782811,0.0611903065302504,78,0.0201380741457537,0.0402741065308831,0.0604060556017787,79,0.0198831702670611,0.0397643753301354,0.059641650179499,80,0.0196346386743029,0.0392673849212618,0.0588963464959308,81,0.019392243371958,0.0387826635436474,0.0581694374862219,82,0.0191557598755964,0.038309762442233,0.0574602505521723,83,0.0189249745184607,0.0378482544820363,0.0567681454875653,84,0.018699683807684,0.037397732861809,0.0560925125517379,85,0.0184796938257785,0.036957809918259,0.0554327706788301,86,0.0182648196740942,0.0365281160138496,0.0547883658119943,87,0.0180548849544306,0.036108298501119,0.054158769352271,88,0.0178497212861548,0.0356980207576873,0.0535434767132445,89,0.0176491678558172,0.0352969612862607,0.0529420059729848,90,0.017453070996745,0.0349048128745662,0.0523538966157477,91,0.0172612837964412,0.0345212818107782,0.0517787083566959,92,0.0170736657294776,0.0341460871499508,0.0512160200430905,93,0.0168900823141926,0.0337789600278608,0.0506654286263748,94,0.0167104047912038,0.0334196430186019,0.050126548199817,95,0.0165345098222458,0.0330678895327165,0.0495990090968606,96,0.0163622792078739,0.0327234632529711,0.0490824570458204,97,0.016193599622553,0.0323861376050035,0.0485765523768222,98,0.0160283623660578,0.0320556952603521,0.0480809692772387,99,0.0158664631298188,0.0317319276696139,0.0475953950922138,100,0.0157078017774206,0.0314146346236388,0.0471195296672208)
dm <- t(matrix(dv,4))
data <- list(n=dm[,1], d1=dm[,2], d2=dm[,3], d3=dm[,4])


# optim works
f<-function(par, x, y){
  y_p <- par[1] * sin(par[2]*pi/x)
  y_m <- y
  return(sum((y_m - y_p)^2))
}
optim(par=c(1,1), f, method='BFGS', x=data$n, y=data$d2, control = list(maxit=10^4))
#note the extension of maximum n of itterations

# nls fails (even using the ideal starting conditions)
nls(d2 ~ b * sin(a * pi/x), start = list(a=0.5, b=2), data=data)
# it is a problem with convergence and you can force it to show output anyway 
nls(d2 ~ b * sin(a * pi/x), start = list(a=0.5, b=2), data=data, control = nls.control(warnOnly=TRUE))
# or it works without that hack when you add a tiny bit of noise
nls(d2+rnorm(98,0,0.00001) ~ b * sin(a * pi/x), start = list(a=0.5, b=2), data=data)


# residuals
plot(data$d1-2*sin(0.25*pi/x) )
plot(data$d2-2*sin(0.5 *pi/x) )
plot(data$d3-2*sin(0.75*pi/x) )
$\endgroup$
3
  • $\begingroup$ Weterings - I think you got it right, because it makes sense with my experiment. "I found this function by realizing that your values resemble much the length of sides of polygons whose vertices are on a unit circle" But what to do if you dont have any idea on the type of function? --- yes my nls try failed as you also mention $\endgroup$ – pma Sep 25 '17 at 16:13
  • $\begingroup$ plotting the residuals with plot(data$d1-2*sin(0.25*pi/dat$n) ) shows that residuals increase with n. May this be a result of rounding errors in the data I provided? so smaller numbers have greater errors? $\endgroup$ – pma Sep 25 '17 at 17:01
  • $\begingroup$ The data gathering method is also (meta)data. It is likely the best source of information regarding the intruiging $\varpropto n$ scaling of the error-terms (and why the third column has much larger errors). So looking at this additional data is what you do when you have no idea about the function. More strictly, there can be no analysis without having some idea about the function. If you don't have any additional information then numbers are just numbers, and not data. Also without well defined requirements we can do anything with the numbers, a spline with 98 knots fits your data perfectly. $\endgroup$ – Sextus Empiricus Sep 27 '17 at 11:08
1
$\begingroup$

The functions look like reciprocal functions, and indeed, when I run

with(thedata, optimize(function(p) sum((1/(p*n) - d1)^2), c(.0001, 10)))

I get close fits with a p of .6399901 for d1, .3251379 for d2, and .2225606 for d3. The fit is not quite as close for d2 and d3 as d1, so as the parameter distinguishing the dNs varies, there's probably some departure from a pure reciprocal function times a constant as in the case of d1. It would be hard to guess exactly what's going on without more dNs. Experimenting with d2 and d3 suggests the exponent is changing from -1.

$\endgroup$
1
$\begingroup$

As a first guess at the functional form of your data, try plotting $\log(y)$ against $\log(x)$. If this graph looks like a straight line, say $$ \log y = m \log x + b\tag1 $$ then the relationship between $y$ and $x$ is: $$ y = K x^m\tag2 $$ for some constant $K$ (which equals $e^b$). In R, you can fit the form (1) to your logged data using lm() and deduce the parameters $K$ and $m$ from the coefficients returned by the linear fit.

For example, applying (1) to your d1 data yields:

dat <- read.csv("mydata.dat", sep=";", header=T)

fit1 <- lm(log(d1) ~ log(n), data=dat)
b <- fit1$coefficients[1] # intercept
m <- fit1$coefficients[2] # slope
K <- exp(b)

# Result:
# b = 0.4468282
# m = -0.9988257 
# K = 1.563346

Alternatively, you can fit the functional form (2) to your original data using the nls() function. In that case you can either treat $m$ as a parameter to be fitted, or assume $m=-1$ based on the fit returned by lm().

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.