{{{id=19| %auto # reset() # commented out for now -- ticket 7255 forget() # special equation rendering def render(x,name = "temp.png",size = "normal"): if(type(x) != type("")): x = latex(x) latex.eval("\\" + size + " $" + x + "$",{},"",name) var('a b c d e f g h j k l x y z') class Pair: def __init__(self,d): self.x = d[0] self.y = d[1] def __repr__(self): return "%e,%e" % (self.x,self.y) def show_mat(data): yt = len(data) xt = len(data[0]) for y in range(yt): for x in range(xt): sys.stdout.write("%12g " % data[y][x]) print("") print("****") def show_mat_latex(data,fn): yt = len(data) xt = len(data[0]) ya = [] for y in range(yt): xa = [] for x in range(xt): xa.append("%12g " % data[y][x]) ya.append(join(xa," & ")) ys = join(ya,"\\\\") # specify right justification cstr = "r" * xt mat = "$\\left( \\begin{array}{%s}" % cstr mat += ys mat += "\\end{array}\\right)$" render(mat,fn) print "" # example data data1 = [ [-1,-1], [0,3], [1,2.5], [2,5], [3,4], [5,2], [7,5], [9,4] ] def set_data(x): return [Pair(n) for n in x] # polynomial regression of order "p" def polyregress(data,p,debug = False): # matrix rows and columns = p+1 p += 1 n = len(data) # precalculate matrix cell data array mda = [n] + [sum(data[a].x^b for a in range(n)) for b in range(1,2*p-1)] # create/populate square matrix with appended column m = [[mda[r+c] for c in range(p)] + [0.0] for r in range(p)] if debug: show_mat_latex(m,"nmat2.png") # populate RH appended column for v in data: m[0][p] += v.y for r in range(1,p): m[r][p] = sum(v.x^r*v.y for v in data) if debug: show_mat_latex(m,"nmat3.png") # reduce matrix mm=Matrix(RR,m) mm.echelonize() if debug: show_mat_latex(mm.rows(),"nmat4.png") # extract result column a = mm.column(p).list() # build term list y = 0 for j in range(p): y += a[j]*x^j return y # correlation coefficient (r^2) def corr_coeff(data,fy): r = 0 sx = sx2 = sy = sy2 = sxy = 0 n = len(data) for v in data: x = fy(x=v.x) y = v.y sx += x sy += y sxy += x * y sx2 += x * x sy2 += y * y divisor = sqrt((sx2 - (sx*sx)/n) * (sy2 - (sy*sy)/n)) if(divisor != 0): r = ((sxy-(sx*sy)/n)/divisor)^2 return r # standard error def std_error(data,fy): r = 0 n = len(data); if(n > 2): a = 0 for v in data: a += (fy(x=v.x) - v.y)^2 r = sqrt(a/(n-2)) return r # data range def list_minmax(data): max = -1e6 min = 1e6 for v in data: if(max < v.x): max = v.x if(min > v.x): min = v.x return(min,max) /// }}} {{{id=81| # interactive version data = set_data(data1) mp = len(data)-1 @interact def _(p = (0..mp)): y = polyregress(data,p) dp = list_plot([[data[n].x,data[n].y] for n in range(len(data))],rgbcolor='red') lp = plot(y,x,list_minmax(data)) cc = corr_coeff(data,y) se = std_error(data,y) lbl = text("correlation coefficient (r^2): %f\nStandard error: %f" % (cc,se),(5,-2),rgbcolor='black',fontsize=10) show(dp+lp+lbl,ymin=-5,ymax=7) ///
p 
}}} {{{id=109| # static version data = set_data(data1) p = 4 y = polyregress(data,p) dp = list_plot([[data[n].x,data[n].y] for n in range(len(data))],rgbcolor='red') lp = plot(y,x,list_minmax(data)) cc = corr_coeff(data,y) se = std_error(data,y) lbl = text("correlation coefficient (r^2): %f\nStandard error: %f" % (cc,se),(5,-2),rgbcolor='black',fontsize=10) show(dp+lp+lbl,ymin=-5,ymax=7,figsize=(5,4)) /// }}} {{{id=96| # generate result matrices var('y') data = set_data(data1) q = polyregress(data,4,True) #print corr_coeff(data,y) #print std_error(data,y) render(y == q,"poly_result1.png") /// }}} {{{id=95| # generate matrix diagram def add_to_mat(a,b): return a + "\\left( \\begin{array}{%s} " % cspec + join(b,"\\\\") + "\\end{array} \\right)" p = 4 ya = [] ltx = "" cspec = "l" * (p+1) for y in range(p+1): xa = [] for x in range(p+1): if(y == 0 and x == 0): xa.append("n") else: ss = "^{" + str(x+y) + "}" if(x+y < 2): ss = "" xa.append("\\sum{x_i}%s" % ss) ya.append(join(xa," & ")) ltx = add_to_mat(ltx,ya) xa = [] for y in range(p+1): xa.append("a_%d" % y) ltx = add_to_mat(ltx,xa) xa = [] xa.append("\\sum{y_i}") for y in range(1,p+1): ss = "^{" + str(y) + "}" if(y < 2): ss = "" xa.append("\\sum{x_i}%s {y_i}" % ss) ltx = add_to_mat(ltx + " = ",xa) # show(html("$" + ltx + "$")) render(ltx,"equ_mat.png") /// }}} {{{id=98| render("$(x_1,y_1),(x_2,y_2),(x_3,y_3), ... ,(x_n,y_n)$","dataset.png","large") /// }}} {{{id=99| render("$y = a_0 + a_1 x + a_2 x^2 + a_3 x^3 ... a_p x^p, p < n$","polyequ.png","large") /// }}} {{{id=100| render("$R_i = y_i - a_0 - a_1 x_i - a_2{x_i}^2 ... -a_p {x_i}^p$","residual.png","large") /// }}} {{{id=101| xvals = [-1,0,1,2,3,5,7,9] for p in range(9): print "%d: %8d" % (p,sum(xvals[n]^p for n in range(len(xvals)))) /// 0: 8 1: 26 2: 170 3: 1232 4: 9686 5: 79256 6: 665510 7: 5686952 8: 49208966 }}} {{{id=102| render("$\\sum{x_i}^{2}$","equ_sumsq1.png") /// }}} {{{id=103| render("$\\sum_{i = 0}^{n-1}{x_i}^{2}$","equ_sumsq2.png") /// }}} {{{id=110| render("$\\sum_{i = 0}^{n-1}{x_i}^{(r+c)}$","equ_sumsq3.png") /// }}} {{{id=107| # test assumption about RH column data = set_data(data1) s = 0 p = 4 for a in range(p+1): s = 0 for n in range(len(data)): s += data[n].x^a * data[n].y show(s) ///
24.5000000000000
106.500000000000
676.500000000000
5032.50000000000
39904.5000000000
}}} {{{id=108| /// }}}