I have a list of observation data {x,y} and its densityplot is as follows. The blue solid line is the linear least squares fitting with unique weighting. Theoretically, x and y should be linearly related. The linear fitting (blue line) on the plot seems not align with the data with maximum counts.

I expect a linear fitting would be something like the red dashed line. I want to perform a linear least sqaures fitting with a weighting inversely proportional to y error over standard deviation, for example,

for 30 <= x < 31, 1/w_i = Abs[y_i - mean[y]]/SD[y]

I tried to make a weight list using the above criteria but it was extremely time-consuming. Is there any built in function for achieving my purpose?

I also upload the data here if anyone would like to take a look. https://drive.google.com/file/d/0B3sHo6MIAjziaTloRnlGbm1Nam8/view?usp=sharing

Many thanks in advance.

enter image description here

  • $\begingroup$ How do you define "y error" and "standard deviation" here? $\endgroup$ – MarcoB Feb 13 '17 at 5:47
  • $\begingroup$ Is your data censored or truncated around -16.1363? (That's the smallest value in the data.) Is that some threshold for the measuring instrument? In other words, why is the bottom so flat? $\endgroup$ – JimB Feb 13 '17 at 7:00
  • $\begingroup$ Yes, the {x, y} data is logarithmic and coded in 8 bits, so -16.1363 is the smallest value, i.e. 10Log[0.0243]. $\endgroup$ – wkong Feb 14 '17 at 2:34

In your previous question on this dataset, you received a great answer from Anton Antonov that used quantile regression to reject outliers in your data via a statistically rigorous method.

Using his QuantileRegression package as he proposed there, the resulting linear fit on the remaining points seems to fit your current request quite well:

raw = Import["your/data/file.txt", "CSV"];
data2 = Map[Reverse, RandomSample[raw, 10000]];
qs = {0.1, 0.9};
qFuncs = QuantileRegression[data2, 20, qs];
pred2 = Map[
           #[[2]] >= qFuncs[[1]][#[[1]]] && #[[2]] <= qFuncs[[2]][#[[1]]] &, 
           Reverse /@ raw

in = Pick[raw, pred2];
out = Complement[raw, in];

inFit = LinearModelFit[in, {1, x}, x];
allFit = LinearModelFit[raw, {1, x}, x];

  {Style[allFit[x], Blue], Style[inFit[x], Red, Dashed]}, {x, 10, 60},
  PlotStyle -> Thickness[0.01], PlotRange -> {-22, 7},
  Prolog -> {Gray, Point@out, Green, Point@in},
  Axes -> False, Frame -> True, ImageSize -> Large

Mathematica graphics

If this is not acceptable, then perhaps you should further clarify your question.

  • $\begingroup$ Thank you so much! Didn't realize that @Anton Antonov 's QuantileRegression package can be used in this way. $\endgroup$ – wkong Feb 14 '17 at 2:28
  • $\begingroup$ When I run the code you suggested above, the evaluation abort and the Kernel quit occasionally. It did not happen every time, maybe around once every 5 trial. When I tried other data sets, this situation happened more frequently in some of the data sets. Do you have any suggestion for me please? Thanks!! $\endgroup$ – wkong Feb 16 '17 at 2:54

This is an extended comment. First, to agree with @MarcoB, Anton Antonov gave you a great answer. Second, despite the theory stating things are linear, your data provides no support for that theory (at least not across the whole range of the predictor variable).

Because your data has multiple observations for each predictor value, it makes sense to plot the means and confidence intervals associated with each unique predictor value. That plot is created below and maybe things are approximately linear from the predictor values from about 33 and higher, the remaining data supports a separate and smaller slope. Also note that the confidence intervals are so small that only for values below 20 and above about 55 can they even be distinquished from the red point representing the mean.

Quantile regression bends with your data. Try it.

plot of data and means


I used R to create the above plot. Here are the Mathematica commands to produce an equivalent plot:

se = StandardDeviation[#]/Sqrt[Length[#]] & /@ GatherBy[data, First];
means = Mean[#] & /@ GatherBy[data, First];
lower = means - 1.96 se;
upper = means + 1.95 se;
Show[ListPlot[{data, means}, 
  PlotStyle -> {Automatic, {PointSize[0.01], Red}}],
 ListPlot[Transpose[{lower, upper}], Joined -> True, PlotStyle -> Red]]

Data and means by predictor value

2nd Update

Because you've now stated that the data is essentially censored so that the dependent variable values can't go below -16.1363, your data are censored. This suggests that the median might be a better descriptor of the relationship in that it is less influenced by the censoring.

medians = Median[#] & /@ GatherBy[data, First];
means = Mean[#] & /@ GatherBy[data, First];
ListPlot[{data, means, medians}, 
 PlotStyle -> {Automatic, {PointSize[0.01], Red}, {PointSize[0.01], 
    Blue}}, PlotLegends -> {"Raw data", "Means", "Medians"}]

Data, means, and medians

  • $\begingroup$ Thank you for your comment. May I know which function you have used to make the plot? Thanks! $\endgroup$ – wkong Feb 14 '17 at 2:35
  • $\begingroup$ My list manipulation skills in Mathematica are not strong so I used R to produce the means and the plot. $\endgroup$ – JimB Feb 14 '17 at 5:42
  • $\begingroup$ Actually, I guess I do know how to get the plot of the means in Mathematica: ListPlot[{data, Mean[#] & /@ GatherBy[data, First]}, PlotStyle -> {Automatic, {PointSize[0.02], Red}}] $\endgroup$ – JimB Feb 14 '17 at 6:06
  • $\begingroup$ @wkong I know that more than two years have passed since the above post but, could you post the original data set by providing some other web link? The "drive.google.com/file/d/0B3sHo6MIAjziaTloRnlGbm1Nam8/…" link that you provided originally is presently rejected by my web browser as presenting some sort of risk. Thank you. $\endgroup$ – Gilmar Rodriguez Pierluissi Mar 13 '19 at 20:47
  • $\begingroup$ @GilmarRodriguezPierluissi I am so sorry. The file has been removed from the my Google Drive a long time ago and the original file is already missing. $\endgroup$ – wkong Mar 29 '19 at 2:03

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.