gnuplot-info

[Gnuplot-info] Help request - linear regression giving screwy results

[Eric R. <es...@th...>]

I've been attempting to do trend analysis on U.S. smartphone marketshare.
You can see some of my successful visualizations at 

    http://www.catb.org/esr/comscore/

The problem is that I tried to do a linear-regression fit on Android's
market share and got screwy results.  

To reproduce, save the enclosed comscore.dat and comscore-report, then
do "comscore-report -p Android foo.png" and display foo.png.  What I
see is that the linear-fit line passes through one datapoint but is nowhere
near close to the right slope.  I don't get an error message.

Thanks in advance,
-- 
		<a href="http://www.catb.org/~esr/">Eric S. Raymond</a>

Re: [Gnuplot-info] Help request - linear regression giving screwy results

[Thomas S. <t.s...@fz...>]

on the x-axis you have seconds since 2000-01-01, so your x-values
are around 315360000.
your y-values are around 20.

this leads to parameters which are different by at least 7 orders of
magnitude.

to bring the parameters into the same range, subtract 10 years and scale the
resulting x-values:

offset=10*365*24*60*60
f(x)=1.e-7*m*(x-offset)+b
m=7.
b=6.
fit f(x) 'comscore.dat' using 1:2 via m,b


Eric Raymond-3 wrote:
> 
> I've been attempting to do trend analysis on U.S. smartphone marketshare.
> You can see some of my successful visualizations at 
> 
>     http://www.catb.org/esr/comscore/
> 
> The problem is that I tried to do a linear-regression fit on Android's
> market share and got screwy results.  
> 
> To reproduce, save the enclosed comscore.dat and comscore-report, then
> do "comscore-report -p Android foo.png" and display foo.png.  What I
> see is that the linear-fit line passes through one datapoint but is
> nowhere
> near close to the right slope.  I don't get an error message.
> 
> Thanks in advance,
> -- 
> 		 http://www.catb.org/~esr/ Eric S. Raymond 
> 
> #!/usr/bin/env python
> """
> Generate reports from raw comScore market-share data
> 
> -s = tabulate or plot marketshare trends
> -u = tabulate or plot userbase trends
> -d = tabulate or plot changes in userbase by month
> -w = generate HTML table to stdout; without this, make a plot to a file
> -t = generate text table to stdout
> 
> The raw data is assumed to be in comscore.dat.
> """
> 
> import os, sys, getopt, tempfile, copy
> 
> class comScore:
>     def __init__(self, data):
>         self.data = data
>     #
>     # Framework code
>     #
>     def arithmetize(self):
>         "Turn data to numeric, excluding top row, left column, and -
> entries."
>         for i in range(1, len(self.data)):
>             for j in range(1, len(self.data[0])):
>                 if self.data[i][j] != '-':
>                     self.data[i][j] = float(self.data[i][j])
>     def unarithmetize(self):
>         "Turn numeric table data back to strings."
>         w = len(self.data[0])
>         d = len(self.data)
>         for i in range(1, d):
>             for j in range(1, w):
>                 if self.data[i][j] != '-':
>                     self.data[i][j] = "%.2f" % self.data[i][j]
>     def emit(self):
>         "Ship transformed self.data to a file for plotting."
>         (h, name) = tempfile.mkstemp()
>         ofp = open(name, "w")
>         d = len(self.data)
>         for i in range(d):
>             ofp.write("\t".join(self.data[i]) + "\n")
>         ofp.close()
>         return name
>     def textize(self, ofp=sys.stdout):
>         "Dump data as a tab-separated-values file."
>         for i in range(len(self.data)):
>             ofp.write("\t".join(self.data[i]) + "\n")
>     def webize(self, ofp=sys.stdout):
>         "Generate a table suitable for web display from specified
> self.data."
>         d = len(self.data)
>         w = len(self.data[0])
>         for i in range(d):
>             self.data[i][0] = self.data[i][0][:3] + " " +
> self.data[i][0][3:]
>         ofp.write("<table border='1'>\n")
>         for j in range(w):
>             ofp.write("<tr>")
>             for i in range(d):
>                 ofp.write("<td>" + self.data[i][j] + "</td>")
>             ofp.write("</tr>\n")
>         ofp.write("</table>\n")
>     def lastmonth(self):
>         "Return last month for which report is valid."
>         return self.data[len(self.data)-1][0]
>     def select(self, platform):
>         "Select out data for a single platform."
>         i = self.data[0].index(platform)
>         for j in range(len(self.data)):
>             self.data[j] = [self.data[j][0], self.data[j][i]]
>         self.data.pop(0)
>     #
>     # Data reduction
>     #
>     def usercount(self):
>         "Multiply market shares by smartphone userbase size (last
> column)."
>         w = len(self.data[0])
>         d = len(self.data)
>         self.data[0][w-1] = "Total"
>         for i in range(1, d):
>             for j in range(1, w-1):
>                 if self.data[i][j] != '-':
>                     self.data[i][j] *= self.data[i][w-1]
>                     self.data[i][j] /= 100.0
>     def deltas(self):
>         "Turn self.data into a differences table."
>         w = len(self.data[0])
>         d = len(self.data)
>         differences = copy.deepcopy(self.data)
>         for i in range(2, d):
>             for j in range(1, w):
>                 if self.data[i][j] == '-' or self.data[i-1][j] == '-':
>                     differences [i][j] = '-'
>                 else:
>                     differences [i][j] = self.data[i][j] -
> self.data[i-1][j]
>         # Remove first row, for which there is no corresponding delta.
>         self.data = differences[:1] + differences[2:]
> 
> coreplot = """
> set terminal png nocrop enhanced
> set output '%(output)s'
> set key outside right top vertical Right noreverse noenhanced autotitles
> nobox
> set datafile missing '-'
> set style data linespoints
> set xtics border in scale 1,0.5 nomirror rotate by -45 offset character 0,
> 0, 0
> set xtics norangelimit
> set xtics ()
> plot '%(input)s' using 2:xtic(1) title columnheader(2), \
>      		    '' using 3:xtic(1) title columnheader(3), \
>      		    '' using 4:xtic(1) title columnheader(4), \
>      		    '' using 5:xtic(1) title columnheader(5), \
>      		    '' using 6:xtic(1) title columnheader(6)
> """
> 
> predictive = """
> set terminal png nocrop enhanced
> set output '%(output)s'
> set datafile missing '-'
> set style data points
> set xtics border in scale 1,0.5 nomirror rotate by -45 offset character 0,
> 0, 0
> set xtics norangelimit
> set xtics ()
> set xdata time
> set timefmt '%%b%%Y'
> set xtics format '%%b%%Y'
> unset key
> f(x) = m*x + b
> fit f(x) '%(input)s' using 1:2 via m,b
> plot '%(input)s' using 1:2, f(x)
> """
> 
> def gnuplot(inputname, plot, outputname):
>     "Generate a derived plot."
>     plot = plot % {"output" : outputname, "input" : inputname}
>     print plot
>     ofp = os.popen("gnuplot -", "w")
>     ofp.write(plot)
>     ofp.close()
> 
> def grab(filename):
>     "Grab the contents of a data file."
>     lines = []
>     for line in open(filename):
>         if line[0] != '#':
>             lines.append(line.strip().split("\t"))
>     return lines
>  
> if __name__ == '__main__':
>     (options, arguments) = getopt.getopt(sys.argv[1:], "suwdtmp:")
>     plotprefix = ""
>     tabulate = False
>     textdump = False
>     datedump = False
>     share = False
>     user = False
>     deltas = False
>     predict = None
> 
>     basedata = grab('comscore.dat')
>     info = comScore(basedata)
> 
>     for (opt, val) in options:
>         if opt == '-s':
>             share = True
>         elif opt == '-u':
>             user = True
>         elif opt == '-d':
>             deltas = True
>         elif opt == '-w':
>             tabulate = True
>         elif opt == '-t':
>             textdump = True
>         elif opt == '-m':
>             datedump = True
>         elif opt == '-p':
>             predict = val
> 
>     if user:
>         title = "Userbase by platform, "
>     else:
>         title = "Market-share per platform, "
>     if deltas:
>         title += "change per month, "
>     if user:
>         yformat = "set format y '%%.0fM'\n"
>         title += "units of 1M users."
>     else:
>         yformat = "set format y '%%.0f%%%%'\n"
>         title += "units of 1%%."
>     title = "set title '%s'\n" % title
> 
>     if user or deltas:
>         info.arithmetize()
>         if user:
>             info.usercount()
>         if deltas:
>             info.deltas()
>         info.unarithmetize()
> 
>     if tabulate:
>         info.webize()
>     elif textdump:
>         info.textize()
>     elif datedump:
>         sys.stdout.write(info.lastmonth())
>     elif predict:
>         info.select(predict)
>         title = 'set title "%s"\n' % predict
>         gnuplot(info.emit(), title + yformat + predictive, arguments[0])
>     else:      
>         gnuplot(info.emit(), title + yformat + coreplot, arguments[0])
> 
> # End
> 
> ------------------------------------------------------------------------------
> Achieve unprecedented app performance and reliability
> What every C/C++ and Fortran developer should know.
> Learn how Intel has extended the reach of its next-generation tools
> to help boost performance applications - inlcuding clusters.
> http://p.sf.net/sfu/intel-dev2devmay
> _______________________________________________
> gnuplot-info mailing list
> gnu...@li...
> https://lists.sourceforge.net/lists/listinfo/gnuplot-info
> 
> 

-- 
View this message in context: http://old.nabble.com/Help-request---linear-regression-giving-screwy-results-tp31627742p31627998.html
Sent from the Gnuplot - User mailing list archive at Nabble.com.

Re: [Gnuplot-info] Help request - linear regression giving screwy results

[Eric R. <es...@th...>]

Thomas Sefzick <t.s...@fz...>:
> 
> on the x-axis you have seconds since 2000-01-01, so your x-values
> are around 315360000.
> your y-values are around 20.
> 
> this leads to parameters which are different by at least 7 orders of
> magnitude.
> 
> to bring the parameters into the same range, subtract 10 years and scale the
> resulting x-values:
> 
> offset=10*365*24*60*60
> f(x)=1.e-7*m*(x-offset)+b
> m=7.
> b=6.
> fit f(x) 'comscore.dat' using 1:2 via m,b

Thank you, that does solve the problem.  I take it the failure of the
fit to converge properly was a result of some internal problem due
to floating-point limitations?

It seems as though this might be a general issue with plots using
time-valued x data. Might I suggest we do a FAQ entry on this topic?
If you were to post a draft here I would read and try to polish it it
for comprehensibility from a user point of view.
-- 
		<a href="http://www.catb.org/~esr/">Eric S. Raymond</a>

Re: [Gnuplot-info] Help request - linear regression giving screwy results

[Thomas S. <t.s...@fz...>]

let's look at your data:

Dec2009 5.2
...
Mar2011 34.7

dates in gnuplot are calculated as 'seconds since beginning
of the gnuplot epoch' which is '2000-01-01 00:00:00'. so internally
your data are:

gnuplot> print strptime("%b%Y","Dec2009")
312940800.0
gnuplot> print strptime("%b%Y","Mar2011")
352252800.0

312940800.0 5.2
...
352252800.0 34.7

this means the slope 'm' of your linear regression function
f(x) = m*x + b
will be around
gnuplot> print (34.7-5.2)/(352252800.0-312940800.0)
7.50407000407001e-07

the y-axis intercept 'b' will be around
gnuplot> print 34.7 - 7.50407000407001e-07*352252800.0
-229.632967032967

thus the fit algorithm should make small steps in the 1.e-8 range
when varying 'm', but large steps in the 1.e-1 range when varying
'b'. this difference in step size makes it (nearly) impossible for the
fit algorithm to converge properly, so it will end up with a regression
function which crosses your data points but with a wrong slope.

a small step (1.e-8) which has a strong effect on 'm' has 'nearly'
no effect on 'b'.
and a step which has an effect on 'b' will - when applied to 'm' -
result in a regression function which misses the data points totally.
this behavior confuses the fit algorithm.

to make the fit algorithm happy, apply a factor of e.g. 1.e-7 to 'm'

f(x) = 1.e-7 * m *x + b

and the fit algorithm will converge after a few iterations because
the effect of a step in 'm' or 'b' on the deviation of the regression
function from the data is now nearly equal.

not necessary, but also helpful: shift the origin of the regression
function into the x-range defined by the data points (=subtract
10 years from 'x'), this makes the fit algorithm faster.


Eric Raymond-3 wrote:
> 
> Thomas Sefzick <t.s...@fz...>:
>> 
>> on the x-axis you have seconds since 2000-01-01, so your x-values
>> are around 315360000.
>> your y-values are around 20.
>> 
>> this leads to parameters which are different by at least 7 orders of
>> magnitude.
>> 
>> to bring the parameters into the same range, subtract 10 years and scale
>> the
>> resulting x-values:
>> 
>> offset=10*365*24*60*60
>> f(x)=1.e-7*m*(x-offset)+b
>> m=7.
>> b=6.
>> fit f(x) 'comscore.dat' using 1:2 via m,b
> 
> Thank you, that does solve the problem.  I take it the failure of the
> fit to converge properly was a result of some internal problem due
> to floating-point limitations?
> 
> It seems as though this might be a general issue with plots using
> time-valued x data. Might I suggest we do a FAQ entry on this topic?
> If you were to post a draft here I would read and try to polish it it
> for comprehensibility from a user point of view.
> -- 
> 		 http://www.catb.org/~esr/ Eric S. Raymond 
> 
> ------------------------------------------------------------------------------
> Achieve unprecedented app performance and reliability
> What every C/C++ and Fortran developer should know.
> Learn how Intel has extended the reach of its next-generation tools
> to help boost performance applications - inlcuding clusters.
> http://p.sf.net/sfu/intel-dev2devmay
> _______________________________________________
> gnuplot-info mailing list
> gnu...@li...
> https://lists.sourceforge.net/lists/listinfo/gnuplot-info
> 
> 

-- 
View this message in context: http://old.nabble.com/Help-request---linear-regression-giving-screwy-results-tp31627742p31635555.html
Sent from the Gnuplot - User mailing list archive at Nabble.com.

Re: [Gnuplot-info] Help request - linear regression giving screwy results

[Eric R. <es...@th...>]

Thomas Sefzick <t.s...@fz...>:
> to make the fit algorithm happy, apply a factor of e.g. 1.e-7 to 'm'
> 
> f(x) = 1.e-7 * m *x + b
> 
> and the fit algorithm will converge after a few iterations because
> the effect of a step in 'm' or 'b' on the deviation of the regression
> function from the data is now nearly equal.
> 
> not necessary, but also helpful: shift the origin of the regression
> function into the x-range defined by the data points (=subtract
> 10 years from 'x'), this makes the fit algorithm faster.

OK, not a floating-point misconvergence but something more subtle.  
Let's see if I can recast this into a FAQ emtry.  Criticize, please.

-----------------------------------------------------------------------------
Q: I'm using time xdata and my curve-fit produces absurd results. How do I 
fix this?

A: Change the coefficients on your fit function so that incrementing any
of them produces a difference in the value of y that is about the same as
incrementing any other.

Say, for example, that you are trying a linear fit on some time-based 
data that looks like this:

2009-12-00 5.2
2011-03-00 34.7

Because dates in gnuplot are calculated as 'seconds since beginning
of the gnuplot epoch' which is '2000-01-01 00:00:00', your data
actually looks to the curve-fitter like this:

312940800.0 5.2
352252800.0 34.7

It's the large difference in magnitude between x and y values that causes 
the problem.  The fit algorithm will make small steps in the 1.e-8 range
when varying 'm', but large steps in the 1.e-1 range when varying
'b'. This difference in step size makes it (nearly) impossible for the
fit algorithm to converge properly, so it will end up with a regression
function which crosses your data points but with a wrong slope.

To fix this, apply a scale factor to the m coefficient that makes values of
mx comparable in magnitude to values of y.  In this case fitting to 
f(x) = 1.e-7 * m *x + b will work much better.

Shifting the origin of the regression function to inside the x-range defined 
by the data points is not necessary for convergence, but it will make the fit 
algorithm run faster.  In this case, try subtracting 10 years, like so:

offset=10*365*24*60*60
f(x)=1.e-7*m*(x-offset)+b

The fix (and the optimization) generalizes to nonlinear fit functions.
The curve-fitting algorithm is going to walk through tuples of
coefficients looking for a tuples with small residuals.  If there are
large enough variations in the magnitude of changes in y as different
coefficience are tweaked, that search may converge on a wrong tuple.
To fix this, tweak the function so the search has roughly equal grain
in all directions.
-----------------------------------------------------------------------------

-- 
		<a href="http://www.catb.org/~esr/">Eric S. Raymond</a>