{VERSION 6 0 "IBM INTEL NT" "6.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 256 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 261 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 266 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 267 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 268 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 271 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 272 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 273 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 } {PSTYLE "Maple Output" 0 11 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Plot" -1 256 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Normal" -1 257 1 {CSTYLE "" -1 -1 "T imes" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 0 0 1 0 1 0 2 2 0 1 } {PSTYLE "Maple Output" -1 258 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 3 0 0 0 0 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 12 "with(stats):" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 106 "Notice that the following commands use a lot o f parenthesis, so be careful when you try your own examples." }}} {EXCHG {PARA 0 "" 0 "" {TEXT 256 9 "Example 1" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "XY1:=[[0,1],[2,2],[4,1]]:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 42 "PlotXY1:=pointplot(XY1):\ndisplay(PlotXY1);" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "fit[leastsquare[[X,Y]]]([[ 0,2,4],[1,2,1]]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 108 "Line1 :=plot(4/3, X=0..4,Y=0..2, color=black):\ndisplay(\{Line1,PlotXY1\}, t itle=`Example 1 least squares fit`);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 257 9 "Example 2" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "XY2:=[[ 0,3],[1,5],[2,5]]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "PlotX Y2:=pointplot(XY2):\ndisplay(PlotXY2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "fit[leastsquare[[X,Y]]]([[0,1,2],[3,5,5]]);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 113 "Line2:=plot(10/3 + X, X=0.. 4,Y=0..6, color=black):\ndisplay(\{Line2,PlotXY2\}, title=`Example 2 l east squares fit`);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 258 9 "Example 3 " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "XY3:=[[0,3],[1,5],[2,6] ,[3,7]]:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "PlotXY3:=pointp lot(XY3):\ndisplay(PlotXY3);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "fit[leastsquare[[X ,Y]]]([[0,1,2,3],[3,5,6,7]]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 123 "Line3:=plot(33/10 + (13/10) *X, X=0..4,Y=0..8, color=black):\nd isplay(\{Line3,PlotXY3\}, title=`Example 3 least squares fit`);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT 261 31 "How do you test for power laws?" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 190 "What if \+ our points don't lie in a straight line? What if it looks like it sho uld be a quadratic function or even a square root function? In other \+ words, we want to see if there is a power " }{TEXT 272 1 "m" }{TEXT -1 32 " and a proportionality constant " }{TEXT 273 1 "b" }{TEXT -1 20 " so that the formula" }}{PARA 257 "" 0 "" {XPPEDIT 18 0 "y = b*x^m ;" "6#/%\"yG*&%\"bG\"\"\")%\"xG%\"mGF'" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 36 "effectively mirrors the real data. " }{TEXT 260 60 "T aking (natural) logarithms of the proposed power law yields" }}{PARA 257 "" 0 "" {XPPEDIT 18 0 "ln(y) = ln(b)+m*ln(x);" "6#/-%#lnG6#%\"yG,& -F%6#%\"bG\"\"\"*&%\"mGF,-F%6#%\"xGF,F," }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 16 "So, if we write " }{XPPEDIT 18 0 "Y = ln(y);" "6#/% \"YG-%#lnG6#%\"yG" }{TEXT -1 7 " and " }{XPPEDIT 18 0 "X = ln(x);" " 6#/%\"XG-%#lnG6#%\"xG" }{TEXT -1 3 ", " }{XPPEDIT 18 0 "B = ln(b);" " 6#/%\"BG-%#lnG6#%\"bG" }{TEXT -1 59 ", this becomes the equation of a \+ line in the new variables " }{XPPEDIT 18 0 "X;" "6#%\"XG" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "Y;" "6#%\"YG" }{TEXT -1 1 ":" }}{PARA 257 "" 0 "" {XPPEDIT 18 0 "Y = mX+B;" "6#/%\"YG,&%#mXG\"\"\"%\"BGF'" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 5 "Thus," }{TEXT 263 1 " " }{TEXT 264 207 "in order for there to be a power law for the original data, t he ln-ln data should (approximately) satisfy the equation of a line, a nd vise verse. If we get a good line fit to the ln-ln data, then the \+ slope " }{XPPEDIT 18 0 "m" "6#%\"mG" }{TEXT 270 75 " of this line is t he power relating the original data, and the exponential " }{XPPEDIT 18 0 "exp(B);" "6#-%$expG6#%\"BG" }{TEXT 265 8 " of the " }{XPPEDIT 18 0 "Y" "6#%\"YG" }{TEXT 271 43 "-intercept is the proportionality co nstant " }{XPPEDIT 18 0 "b" "6#%\"bG" }{TEXT 268 26 " in the original \+ relation " }{XPPEDIT 18 0 "y = b*x^m" "6#/%\"yG*&%\"bG\"\"\")%\"xG%\"m GF'" }{TEXT 269 3 ". " }{TEXT -1 161 "With real data it is not too ha rd to see if the ln-ln data is well approximated by a line, in which c ase the original data is well-approximated by a power law. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 262 45 "Number of genes vs number of cells example: " }{TEXT -1 162 "I used \+ colons after most of these commands to suppress the output. If you wa nt to go back and see what each command is doing, replace the colons w ith semicolons." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 259 0 "" } {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "with(linalg):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 382 "points:=\n[[600000,250],[200000,60],[60000,25 ],[10000,12],[2500,5]]:\n#These are our data points. It says Humans h ave 600,000 genes and 250 cell types, Annelid worms have 200,000 genes and 60 cell types, Jellyfish have 60,000 genes and 25 cell types, spo nges have 10,000 genes and 12 cell types, and yeast have 2500 genes an d 5 cell types. [Number of genes, number of cell types]." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "datapoints:=pointplot(points ):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "display(datapoints,\n title=`plot of [gene,cell]`);\n" }}}{PARA 0 "" 0 "" {TEXT -1 204 " And now for the ln-ln data. First we put our points into a matrix. T hen we take the natural log of those values. To keep Maple happy, we \+ then get the decimal or floating point values of our log values." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 26 "matrixdata:=evalm(points):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 248 "lndata:=map(ln,matrixdata): #take ln of the ge ne/cell values\ndecimallndata:=map(evalf,lndata): #Get decimal (float ing point) values.\n #This speeds up computations later in the\n #le ast squares fit - otherwise Maple tries working\n #symbolically." }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "lnlnplot:=pointplot(decimall ndata):\ndisplay(lnlnplot,title=`ln-ln data`);\n" }}}{PARA 0 "" 0 "" {TEXT -1 135 "Notice the ln-ln plot really does seem close to a line! \+ How do we find the best line to fit this data, other than just eye-ba lling it?" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 267 0 "" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 266 58 "How do I find the \+ best line fit to a collection of points?" }}{PARA 0 "" 0 "" {TEXT -1 249 " Well, any good calculator or mathematical software already h as a canned program to solve this \"linear regression\", or \"least sq uares\" problem. But it's not hard to explain what's going on: Label our data points for which we seek a line fit as" }}{PARA 257 "" 0 "" {XPPEDIT 18 0 "[[X[1], Y[1]], [X[2], Y[2]], [X[3], Y[3]], `...`, [X[n] , Y[n]]];" "6#7'7$&%\"XG6#\"\"\"&%\"YG6#F(7$&F&6#\"\"#&F*6#F/7$&F&6#\" \"$&F*6#F5%$...G7$&F&6#%\"nG&F*6#F<" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 40 "The solution to our problem is the pair " }{XPPEDIT 18 0 "[m, B];" "6#7$%\"mG%\"BG" }{TEXT -1 94 " for which the the sum of the squared vertical distances between the data points and the line " } {XPPEDIT 18 0 "Y = mX+B" "6#/%\"YG,&%#mXG\"\"\"%\"BGF'" }{TEXT -1 35 " is minimized. Precisely, minimize" }}{PARA 257 "" 0 "" {TEXT -1 1 " \+ " }{XPPEDIT 18 0 "F(m,B) := Sum((m*X[i]+B-Y[i])^2,i = 1 .. n);" "6#>-% \"FG6$%\"mG%\"BG-%$SumG6$*$,(*&F'\"\"\"&%\"XG6#%\"iGF/F/F(F/&%\"YG6#F3 !\"\"\"\"#/F3;F/%\"nG" }}{PARA 0 "" 0 "" {TEXT -1 57 "by setting its d erivatives with respect to the variables " }{XPPEDIT 18 0 "m, B" "6$% \"mG%\"BG" }{TEXT -1 62 " equal to zero. This leads to a system of tw o equations for " }{XPPEDIT 18 0 "m" "6#%\"mG" }{TEXT -1 5 " and " } {XPPEDIT 18 0 "B" "6#%\"BG" }{TEXT -1 2 ", " }}{PARA 257 "" 0 "" {XPPEDIT 18 0 "matrix([[Sum(X[i]^2,i = 1 .. n), Sum(X[i],i = 1 .. n)], [Sum(X[i],i = 1 .. n), n]])*matrix([[m], [B]]) = matrix([[Sum(X[i]*Y[ i],i = 1 .. n)], [Sum(Y[i],i = 1 .. n)]]);" "6#/*&-%'matrixG6#7$7$-%$S umG6$*$&%\"XG6#%\"iG\"\"#/F1;\"\"\"%\"nG-F+6$&F/6#F1/F1;F5F67$-F+6$&F/ 6#F1/F1;F5F6F6F5-F&6#7$7#%\"mG7#%\"BGF5-F&6#7$7#-F+6$*&&F/6#F1F5&%\"YG 6#F1F5/F1;F5F67#-F+6$&FU6#F1/F1;F5F6" }{TEXT -1 1 "." }}{PARA 258 "" 1 "" {TEXT -1 109 "This system always has a unique solution, and your \+ calculator solves it to find the linear regression line. " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "with(stats): #first load the statis tics commands. " }}}{PARA 0 "" 0 "" {TEXT -1 95 "Now we need to pick o ff the x-values and the y-values to get them into the least-square syn tax." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 295 "Xs:=convert(col(dec imallndata,1),list): #convert the first column of\n #the ln-ln data i nto a list of the \"x's\" The MAPLE least squares\n #command wants t o have lists as the input, not matrix columns, even\n #though it's ha rd for us to see any difference\nYs:=convert(col(decimallndata,2),list ):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "fit[leastsquare[[X,Y] ]]([Xs,Ys]);#this is the bizarre syntax" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 256 "" 0 "" {TEXT -1 65 "We can paste in th e equation of the line and see how well we did." }}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 125 "line:=plot(-3.721684795+.6636653414*X, X=5..1 5,Y=0..6,\n color=black):\ndisplay(\{line,lnlnplot\}, title=`least s quares fit`);\n" }}}{PARA 0 "" 0 "" {TEXT -1 70 "Finally, we can go ba ck from the least squares line fit to a power law" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 72 "m:=.6636653414: #power\nb:=exp(-3.721684795): #proportionality constant\n" }}}{PARA 0 "" 0 "" {TEXT -1 40 "Notice \+ in the command lines below how to" }}{PARA 0 "" 0 "" {TEXT -1 20 " \+ (i) make a title" }}{PARA 0 "" 0 "" {TEXT -1 81 " (ii) get the axes labeled \"g\" and \"c\" for number of genes and number of cells" }} {PARA 0 "" 0 "" {TEXT -1 210 " (iii) get the display to include app ropriate ranges - to contain all our data\n (iv) first create a pic ture of the power function, and then display its graph along with the \+ original gene/cell data points." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "powerplot:=plot(b*g^m,g=0..700000,c=0..300,color=black):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 174 "display(\{powerplot,datapoi nts\},title=\n `power law for gene-cell data`);\n#by calling the vari ables g and c, and giving their ranges, I\n#get Maple to label the axe s as I want\n" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{PARA 11 "" 1 "" {TEXT -1 0 "" }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}}{MARK "0 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 } {PAGENUMBERS 0 1 2 33 1 1 }