2
$\begingroup$

I have been trying to find an expression that properly models the following plot in Mathematica v9.0:

enter image description here

I have tried modelling it as odd-degree polynomials using Fit, though they are not as true to the plot as I would like. I have also tried some rational functions using FindFit, though they did not produce great results either. Here is the data for the plot (please excuse the length):

data = {{0.01, 31323.911983664115}, {0.02, 7828.515278216496}, 
{0.03, 3477.466902377181}, {0.04, 1954.6003707196044}, 
{0.05, 1249.805747045021}, {0.060000000000000005, 867.0851482013874}, 
{0.06999999999999999, 636.4673453947054}, {0.08, 486.9311903820745}, 
{0.09, 384.53291484195296}, {0.09999999999999999, 311.3866049290406}, 
{0.11, 257.3412596538772}, {0.12, 216.28889448926301}, 
{0.13, 184.3763637628718}, {0.14, 159.07620675640447}, 
{0.15000000000000002, 138.6752691568414}, {0.16, 121.97942143701592}, 
{0.17, 108.13598033671094}, {0.18000000000000002, 96.52312083892006}, 
{0.19, 86.67888692047772}, {0.2, 78.25439415303052}, 
{0.21000000000000002, 70.98224622462985}, {0.22, 64.65476812902567}, 
{0.23, 59.108720309872375}, {0.24000000000000002, 54.21438051092778}, 
{0.25, 49.86762440555636}, {0.26, 45.98410019837538}, 
{0.27, 42.494888103186845}, {0.28, 39.343227748397005}, 
{0.29000000000000004, 36.482023675911016}, {0.3, 33.87192458958273}, 
{0.31, 31.479830384860904}, {0.32, 29.277721415804766}, 
{0.33, 27.24173281624233}, {0.34, 25.351416831221186}, 
{0.35000000000000003, 23.589150578565473}, 
{0.36000000000000004, 21.939657157836862}, {0.37, 20.389615718192935}, 
{0.38, 18.9273417881266}, {0.39, 17.542523416230985}, 
{0.4, 16.226001865111854}, {0.41000000000000003, 14.969588018994479}, 
{0.42000000000000004, 13.76590750924032}, {0.43, 12.608268975384776}, 
{0.44, 11.490550967730906}, {0.45, 10.407103838241714}, 
{0.46, 9.352663616470208}, {0.47000000000000003, 8.322275368773694}, 
{0.48000000000000004, 7.3112239233869944}, {0.49, 6.31497013423562}, 
{0.5, 5.329091069609141}, {0.51, 4.349222660123271}, 
{0.52, 3.371003432129057}, {0.53, 2.3900179927475245}, 
{0.54, 1.4017389237218636}, {0.55, 0.401465682373002}, 
{0.56, -0.6157410034724049}, {0.5700000000000001, 
-1.6551309297530246}, {0.5800000000000001, -2.7223443206129017}, 
{0.59, -3.8235020827272126}, {0.6, -4.965309067111152}, 
{0.61, -6.155173537655297}, {0.62, -7.40134674415246}, 
{0.63, -8.713087397402443}, {0.64, -10.100857004827224}, 
{0.65, -11.576553527513198}, {0.66, -13.153792773163353}, 
{0.67, -14.848249493672865}, {0.68, -16.678073516969427}, 
{0.6900000000000001, -18.664400696558992}, 
{0.7000000000000001, -20.8319844100309}, {0.7100000000000001, 
-23.209981347431455}, {0.72, -25.832936214168225}, 
{0.73, -28.742024906000147}, {0.74, -31.986636416120277}, 
{0.75, -35.62640275803951}, {0.76, -39.73382737038074}, 
{0.77, -44.39772166886549}, {0.78, -49.72774571490505}, 
{0.79, -55.86047672753066}, {0.8, -62.96762146313412}, 
{0.81, -71.26728325169432}, {0.8200000000000001, -81.03965543720359}, 
{0.8300000000000001, -92.64924990010493}, 
{0.8400000000000001, -106.57697657528723}, {0.85, -123.4674219935693}, 
{0.86, -144.20020085637648}, {0.87, -170.00058481852963}, 
{0.88, -202.61642351772565}, {0.89, -244.61139750885124}, 
{0.9, -299.871845000486}, {0.91, -374.5269899793762}, 
{0.92, -478.7211890414361}, {0.93, -630.2802195655562}, 
{0.9400000000000001, -862.9998749185709}, 
{0.9500000000000001, -1247.6385388995275}, 
{0.9600000000000001, -1953.8265188616788}, 
{0.97, -3477.3449030472225}, {0.98, -7828.512945138938}, 
{0.99, -31323.91198390073}}

Any help is very much appreciated.

Improve this question
$\endgroup$
7
  • $\begingroup$ Have you tried FindFormula? $\endgroup$ – Carl Lange May 24 '19 at 19:11
  • $\begingroup$ @CarlLange Unfortunately, I'm working in Mathematica 9.0, which does not support that function. I've edited the question to include that information. $\endgroup$ – Spencer Keller May 24 '19 at 19:12
  • 2
    $\begingroup$ This looks suspiciously like Cot... $\endgroup$ – Carl Lange May 24 '19 at 19:35
  • $\begingroup$ @CarlLange Thanks for the suggestion. I promise I'm not trolling! $\endgroup$ – Spencer Keller May 24 '19 at 19:52
  • $\begingroup$ Oh, don't worry! I didn't mean to insinuate that at all. I wasn't able to get a satisfactory fit with Cot, but it really does seem very close in shape. Hopefully someone with more maths knowledge than me will be able to help. $\endgroup$ – Carl Lange May 24 '19 at 20:20
4
$\begingroup$

Here's a simple start that does not aim at the ultimate fit but rather at getting the divergences at $x=0$ and $x=1$ right. For better performance, use LinearModelFit instead of Fit, use more fitting terms, and add the measured variances of the data points as weights for the fit.

The main idea is that by scaling the data a bit, as

scaleddata = {#[[1]], #[[1]]^2 (1 - #[[1]])^2 #[[2]]} & /@ data;

the pattern becomes quite clear (see first plot below).

Here's a pretty good fit: the first term $\frac{1}{x^2}-\frac{1}{(1-x)^2}$ takes care of the data divergence at the range ends (you can split it into two separate terms $\frac{1}{x^2}$ and $\frac{1}{(1-x)^2}$ if you like; but this does not help here, as the divergence is symmetric), and the remaining terms $1$, $x$, etc. fine-tune the fit (you can add more terms if you like, such as $\frac{1}{x}$, $\frac{1}{1-x}$, $x^2$, $x^3$ etc. as well as others):

f[x_] = Fit[data, {1/x^2 - 1/(1 - x)^2, 1, x}, x];

2.16276 + 3.13315 (-(1/(1 - x)^2) + 1/x^2) + 6.8812 x

First we plot the scaled data together with a scaled version of the fit, in order to show that we've covered the data divergences at $x=0$ and $x=1$:

P1 = ListPlot[scaleddata];
P2 = Plot[x^2 (1 - x)^2 f[x], {x, 0, 1}];
Show[P1, P2]

enter image description here

And here's an unscaled plot of the data and the fit: hard to see anything because of the divergences,

P1 = ListPlot[data, PlotRange -> All];
P2 = Plot[f[x], {x, 0, 1}, PlotRange -> All];
Show[P1, P2]

enter image description here

The fit residuals show more structure, so we should probably add more terms to the fit:

ListPlot[{#[[1]], #[[2]] - f[#[[1]]]} & /@ data, PlotRange -> All]

enter image description here

Improve this answer
$\endgroup$
7
  • 1
    $\begingroup$ @AntonAntonov the function f[x] is a fit to the original data. Only the plot is scaled, to compress the vertical scale. I'll add more explanations. $\endgroup$ – Roman May 25 '19 at 18:44
  • $\begingroup$ "... the function f[x] is a fit to the original data." -- Ok, my mistake. Please update your plots to use PlotRange->All. (See this image.) $\endgroup$ – Anton Antonov May 25 '19 at 19:22
  • $\begingroup$ This method allowed me to produce exactly the results I was looking for. Thank you for making the answer so approachable and easy to read. $\endgroup$ – Spencer Keller May 27 '19 at 13:24
  • $\begingroup$ @SpencerKeller The relative errors of the fit in this answer are not that good in the middle. (As shown here and here.) But may be that is not that important to you, or it is less important than (1) the fit function being succinct and (2) that the fit definitely can be obtained in V 9.0. $\endgroup$ – Anton Antonov May 27 '19 at 17:02
  • 1
    $\begingroup$ @AntonAntonov to see if the residuals are problematic it would be helpful to have estimated error bars for the measurements. In their absence it's not really possible to say if the residuals are truly big or small; we need a scale for such statements. $\endgroup$ – Roman May 29 '19 at 16:06
1
$\begingroup$

(This probably works with Version 9 answer, which OP wants, but I have not checked...)

There are some significant outliers in question's data. That can be seen using the option setting PlotRange->All:

ListPlot[data, PlotRange->All]

enter image description here

After that observation in many ways this question can be seen as a possible duplicate of "Model for log-inear data". Using the related answer though provided only a passable fit, not a (very) good one.

After some experimentation I found a good set of knots for Quantile Regression that provides a very good fit. The code follows.

Update 2019-05-27: Included the fit from the other answer.

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/MonadicProgramming/MonadicQuantileRegression.m"] 

knots = Sort[
   Union[Range[0, 0.2, 0.02], Range[0.2, 0.8, 0.1], 
    Range[0.8, 1, 0.02]]];
knots = Sort[Union[knots[[2 ;; -2]], {0.01, 0.99}]];
Length[knots]

(* 27 *)

qrObj =
  QRMonUnit[data]⟹
   QRMonEchoDataSummary⟹
   QRMonQuantileRegression[knots, 0.5]⟹
   QRMonFit[{1/x^2 - 1/(1 - x)^2, 1, x}]⟹
   QRMonPlot⟹       
   QRMonErrorPlots["RelativeErrors" -> False]⟹
   QRMonErrorPlots["RelativeErrors" -> True];

enter image description here

The fit function obtained might need some tweaking when used in the domain [0,1] not "just" in [0.01,0.99].

qFunc = (qrObj⟹QRMonTakeRegressionFunctions)[0.5];
PiecewiseExpand[qFunc[x]]

enter image description here

Improve this answer
$\endgroup$
2
  • $\begingroup$ Anton, that's not a fit, that's a cubic interpolation. $\endgroup$ – Roman May 25 '19 at 20:25
  • 1
    $\begingroup$ @Roman It is definitely a fit using Quantile Regression with B-splines and the median regression quantile (in this case.) $\endgroup$ – Anton Antonov May 25 '19 at 22:53

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.