I'm doing three weighted linear regressions, and the second two are giving me substantially the same answer as I calculate offline, with the first being very wrong.

What have I done wrong? As far as I can tell, everything is the same, and just copy/pasted, with different data. I've reproduced everything below, with the black line from PGFPlots, and the red line (desired behaviour) from my offline calculations.

The column U is standard deviation, V is variance.

The three plots

\documentclass{article}

\usepackage{pgfplots} 
\usepackage{pgfplotstable} 
\pgfplotsset{compat=1.13}


\begin{document}

%%%%%%%%%%%%%%%%%%%%%%
%This one doesn't work
%%%%%%%%%%%%%%%%%%%%%%
\begin{tikzpicture}[baseline]
\begin{axis}[
    xlabel={$x$},
    ylabel={$y$},
]
\addplot[   
    only marks,
    error bars/.cd, y dir=both, y explicit,
] 
table[
    x = X,
    y = Y,
    y error = U,
]
{
X   Y   U   V
10  34.893  0.6 0.360
20  31.263  1.134   1.286
30  49.149  1.807   3.265
40  53.969  2.717   7.382
50  68.956  2.228   4.964
}; 
\addplot[]
table[
    x = X,
    y={create col/linear regression={y=Y, variance=V}},
]
{
X   Y   U   V
10  34.893  0.6 0.360
20  31.263  1.134   1.286
30  49.149  1.807   3.265
40  53.969  2.717   7.382
50  68.956  2.228   4.964
};
\addplot[red]
table[
    x = X,
    y = yhat,
]
{
X   yhat
10  33.111
20  40.209
30  47.307
40  54.405
50  61.503
};
\end{axis}
\end{tikzpicture}


%%%%%%%%%%%%%%%%%%%%%%
%This one does work
%%%%%%%%%%%%%%%%%%%%%%    
\begin{tikzpicture}[baseline]
\begin{axis}[
    xlabel={$x$},
    ylabel={$y$},
]
\addplot[   
    only marks,
    error bars/.cd, y dir=both, y explicit,
] 
table[
    x = X,
    y = Y,
    y error = U,
]
{
X   Y   U   V
10  10.147  0.848   0.719
20  15.049  0.768   0.590
30  28.045  1.008   1.016
40  34.733  1.778   3.161
50  42.145  2.386   5.693
}; 
\addplot[]
table[
    x = X,
    y={create col/linear regression={y=Y, variance=V}},
]
{
X   Y   U   V
10  10.147  0.848   0.719
20  15.049  0.768   0.590
30  28.045  1.008   1.016
40  34.733  1.778   3.161
50  42.145  2.386   5.693
};
\addplot[red]
table[
    x = X,
    y = yhat,
]
{
X   yhat
10  8.813
20  17.356
30  25.899
40  34.441
50  42.984
};
\end{axis}
\end{tikzpicture}


%%%%%%%%%%%%%%%%%%%%%%
%This one does work
%%%%%%%%%%%%%%%%%%%%%%    
\begin{tikzpicture}[baseline]
\begin{axis}[
    xlabel={$x$},
    ylabel={$y$},
]
\addplot[   
    only marks,
    error bars/.cd, y dir=both, y explicit,
] 
table[
    x = X,
    y = Y,
    y error = U,
]
{
X   Y   U   V
10  16.983  0.417   0.174
20  24.753  0.977   0.955
30  39.675  1.796   3.226
40  38.524  3.838   14.730
50  48.147  3.02    9.120
}; 
\addplot[]
table[
    x = X,
    y={create col/linear regression={y=Y, variance=V}},
]
{
X   Y   U   V
10  16.983  0.417   0.174
20  24.753  0.977   0.955
30  39.675  1.796   3.226
40  38.524  3.838   14.730
50  48.147  3.02    9.120
};
\addplot[red]
table[
    x = X,
    y = yhat,
]
{
X   yhat
10  16.994
20  25.642
30  34.290
40  42.938
50  51.587
};
\end{axis}
\end{tikzpicture}

\end{document}

The data was analysed offline as follows, using the first data set add an example:

enter image description here

where T is the transpose, and -1 is the inverse.

The weights are simply the inverse square of the standard deviation. The intercept is 26.01, and the slope 0.710; this is the equation of the red line in my first plot.

  • @marmot image added – masher Jul 22 at 12:24
  • @marmot the desired output is the red line. In the other two cases, the actual and desired are close enough. – masher Jul 22 at 13:05
  • @marmot Mentioned. – masher Jul 22 at 13:32

I took the first three columns of the first data set and fit them off-line (using Igor Pro). I did an unweighted fit and a weighted fit using the standard deviations as weightings. Below are the results (black is unweighted and red is weighted).

A plot of data from first set generated using Igor Pro. Black is unweighted and red is weighted.

Neither of the lines exactly agree with either of your results. The weighted line I get is close to what you have from the off-line weighted fit. I suggest that one issue is at play in this case. Using variances versus standard uncertainties as weighting factors exposes the algorithm to round-off errors. By reference to your comment of close enough for the other two plots, did you use standard uncertainties as weightings during the off-line fitting and use variances with PFGPlots? This does matter!

Otherwise, I cannot imagine why the fit fails. I compiled your example on my system and got the same result as you did. I changed the order of the graphs (putting the bad one last) and still had an incorrect fit. Finally, I increased the significance of the variances (e.g. 4.96398) and changed the order of the points. Still I had failure. My only other finding is this: When you change the variance of the last data point to 2, the fit line parallels the external fit line.

\documentclass{article}

\usepackage{pgfplots} 
\usepackage{pgfplotstable} 
\pgfplotsset{compat=1.13}

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%
%This one doesn't work
%%%%%%%%%%%%%%%%%%%%%%
\begin{tikzpicture}[baseline]
\begin{axis}[
    xlabel={$x$},
    ylabel={$y$},
]
\addplot[   
    only marks,
    error bars/.cd, y dir=both, y explicit,
] 
table[
    x = X,
    y = Y,
    y error = U,
]
{
X   Y   U   V
10  34.893  0.6 0.360
20  31.263  1.134   1.286
30  49.149  1.807   3.265
40  53.969  2.717   7.382
50  68.956  2.228   4.964
}; 
\addplot[]
table[
    x = X,
    y={create col/linear regression={y=Y, variance=V}},
]
{
X   Y   U   V
10  34.893  0.6 0.360
20  31.263  1.134   1.286
30  49.149  1.807   3.265
40  53.969  2.717   7.382
50  68.956  2.228   2
};
\addplot[red]
table[
    x = X,
    y = yhat,
]
{
X   yhat
10  33.111
20  40.209
30  47.307
40  54.405
50  61.503
};
\end{axis}
\end{tikzpicture}

\end{document}

Perhaps you have discovered a bug in PFGPlots?

  • Thanks. I've added how I did my analysis. The weightings are the inverse square of the standard deviation (inverse of the variance). – masher Jul 24 at 9:12
  • @ChristianFeuersänger could this be a bug? – masher Jul 24 at 13:01
  • @masher - In my analysis, the weights are the inverse of the standard uncertainties. The fit is done using a Levenberg-Marquardt chi-squared optimization. The difference between close enough and exactly the same is to be appreciated here. The reason your answers and mine are NOT exactly the same must be due to round-off errors or differences in the underlying algorithm. – Jeffrey J Weimer Jul 24 at 15:24
  • The difference is probably due to you using different weightings to me, not to round-off errors. – masher Jul 25 at 1:10
  • The deviations are squared to become variances. – Jeffrey J Weimer Jul 25 at 2:26

Your Answer

 

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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