How to find fit for this plot?

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

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.

• Have you tried FindFormula? – Carl Lange May 24 '19 at 19:11
• @CarlLange Unfortunately, I'm working in Mathematica 9.0, which does not support that function. I've edited the question to include that information. – Spencer Keller May 24 '19 at 19:12
• This looks suspiciously like Cot... – Carl Lange May 24 '19 at 19:35
• @CarlLange Thanks for the suggestion. I promise I'm not trolling! – Spencer Keller May 24 '19 at 19:52
• 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. – Carl Lange May 24 '19 at 20:20

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]


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]


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


• @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. – Roman May 25 '19 at 18:44
• "... 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.) – Anton Antonov May 25 '19 at 19:22
• 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. – Spencer Keller May 27 '19 at 13:24
• @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. – Anton Antonov May 27 '19 at 17:02
• @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. – Roman May 29 '19 at 16:06

(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]


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];


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]]


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