Linear regression - trend line with pgfplots 2

Question

Happy new year.

I have fitted a trend line of a population-year data with pgfplots. But I get wrong values of slope m and intercept b. By LibreOffice Calc, I got the values: f(x) = 4,478.87 x – 8,966,996.64 (and also by the regression online tool http://www.alcula.com/calculators/statistics/linear-regression/).

\documentclass{standalone}
\usepackage{pgfplots, pgfplotstable}
\pgfplotsset{compat=newest}
\pgfplotstableread{
X Y
2007 19348
2008 23457
2009 32624
2010 34270
2011 46888
2012 45224
2013 53556
2014 55007
2015 49664
}\datatable

\begin{document}
\begin{tikzpicture}
\begin{axis}[legend pos=outer north east]
\addplot [only marks, mark = *] table {\datatable};
\addplot [thick, red] table[
y={create col/linear regression={y=Y}}
] % compute a linear regression from the input table
{\datatable};
\addlegendentry{$y(x)$}
\addlegendentry{%
$\pgfmathprintnumber{\pgfplotstableregressiona} \; x
\pgfmathprintnumber[print sign]{\pgfplotstableregressionb}$}
\end{axis}
\end{tikzpicture}
\end{document}

If I change the X range of 2007-2015 to 1-9 (i.e., if I move the plot to the left), the slope and intercept results are OK.

\documentclass{standalone}
\usepackage{pgfplots, pgfplotstable}
\pgfplotsset{compat=newest}
\pgfplotstableread{
X Y
1 19348
2 23457
3 32624
4 34270
5 46888
6 45224
7 53556
8 55007
9 49664
}\datatable

\begin{document}
\begin{tikzpicture}
\begin{axis}[legend pos=outer north east]
\addplot [only marks, mark = *] table {\datatable};
\addplot [thick, red] table[
y={create col/linear regression={y=Y}}
] % compute a linear regression from the input table
{\datatable};
\addlegendentry{$y(x)$}
\addlegendentry{%
$\pgfmathprintnumber{\pgfplotstableregressiona} \; x
\pgfmathprintnumber[print sign]{\pgfplotstableregressionb}$}
\end{axis}
\end{tikzpicture}
\end{document}

How could I to fix this issue?

Thanks for your help!

Answer 1

I think this is pushing the limits of the PGF math parser precision. You could instead use gnuplot as a backend to do the fitting:

\documentclass[border=5mm]{standalone}
\usepackage{pgfplots, pgfplotstable}
\pgfplotsset{compat=newest}

\usepackage{filecontents}

\begin{filecontents*}{data.dat}
X Y
2007 19348
2008 23457
2009 32624
2010 34270
2011 46888
2012 45224
2013 53556
2014 55007
2015 49664
\end{filecontents*}

\begin{document}
\begin{tikzpicture}
\begin{axis}[
    legend pos=outer north east,
    xticklabel style={/pgf/number format/1000 sep={}},
    legend cell align=left
]
\addplot [only marks, mark = *] table {data.dat};
\addlegendentry{Data}

\addplot +[raw gnuplot, thick, red, mark=none, smooth] gnuplot {
    FIT_LIMIT=1.e-14;
    f(x)=a*x+b;
    fit f(x) 'data.dat' using 1:2 via a,b;
    % Next, plot the function using the x positions from the table
    plot 'data.dat' using 1:(f($1))
    set print "parameters.dat"; % Open a file to save the parameters into$
    print a,b;
       };     
\addlegendentry[]{\pgfplotstableread{parameters.dat}\parameters% Open the file Gnuplot wrote
    \pgfplotstablegetelem{0}{0}\of\parameters \edef\paramA{\pgfplotsretval}% Get first element, save into \paramA
    \pgfplotstablegetelem{0}{1}\of\parameters \edef\paramB{\pgfplotsretval}% Get second element, save into \paramB
     $\pgfmathprintnumber{\paramA} x \pgfmathprintnumber{\paramB}$    
}  
\end{axis}
\end{tikzpicture}
\end{document}