home *** CD-ROM | disk | FTP | other *** search
/ Class of 2001 / ClassOf2001.iso / Scinotebook / scinoteb / matfuncs < prev    next >
Encoding:
Text File  |  1997-05-15  |  17.7 KB  |  657 lines
  1. (* MATRICES *)
  2.  
  3. tciadj2[m_,i_,j_] := 
  4.   Det[Delete[Transpose[Delete[m,i] ],j]]*(-1)^(i+j)
  5.  
  6. tciadjoint[m_] := 
  7.   If[ !(MatrixQ[m] && Dimensions[m][[1]]==Dimensions[m][[2]]), 
  8.     Return["ERROR: selection must be a square matrix"],
  9.     Module[{ ggg },
  10.       ggg[i_,j_]:= tciadj2[m,i,j];
  11.       Transpose[Array[ggg, Dimensions[m]]] 
  12.     ]
  13.   ]
  14.  
  15. tcicat[a_,b_] :=
  16.   Transpose[Join[Transpose[a],Transpose[b]]]
  17.  
  18. tcipermanent[m_] :=
  19.   If[ !(MatrixQ[m] && Dimensions[m][[1]]==Dimensions[m][[2]]), 
  20.     Return["ERROR: selection must be a square matrix"],
  21.     Module[{ perm=0, n=Dimensions[m][[1]], mm=m, pp, jj, kk, ppp },
  22.       pp = Permutations[Table[i,{i,n-1}]];
  23.       For[jj=1, jj <= Length[pp], jj++,
  24.         ppp = pp[[jj]];
  25.         For[kk=1, kk<n, kk++, mm[[kk]] = RotateRight[m[[kk]],ppp[[kk]] ] ];
  26.         perm += Apply[Plus,Apply[Times,mm]];
  27.       ];
  28.       Return[perm]
  29.     ]
  30.   ]
  31.  
  32. tceigen[m_] :=
  33.   Transpose[Eigenvectors[m]]
  34.    
  35. tcitrace[m_] :=
  36.   Sum[m[[i,i]],{i,Length[m]}]
  37.    
  38. tcicharpoly[m_,v_:0] :=
  39.   If[v===0,
  40.     If[MemberQ[Variables[m],X],
  41.       Return["needvars"],
  42.       Det[X*IdentityMatrix[Length[m]]-m]
  43.     ],
  44.     Det[v*IdentityMatrix[Length[m]]-m]
  45.   ]
  46.    
  47. tcicompanion[p_,v_:0] :=
  48.   Module[{tt, vv, len, ss, aa, bb, var},
  49.     If[v===0,
  50.       vv = Variables[p]; If[Length[vv]==1,var=vv[[1]],Return["needvars"]],
  51.       var=v];
  52.     tt = CoefficientList[p,var];
  53.     len = Length[tt];
  54.     ss = tt/tt[[len]];
  55.     ss = Delete[ss,-1];
  56.     aa = IdentityMatrix[len-1];
  57.     bb = Join[Delete[aa,1],{-ss}]; 
  58.     Transpose[bb]
  59.   ]
  60.    
  61. tcihtranspose[m_] :=
  62.   ComplexExpand[Conjugate[Transpose[m]]]
  63.    
  64. tcinorm[a_,t_:2] := 
  65.   Which[
  66.     VectorQ[a], Which[
  67.       NumberQ[t],     Apply[Plus,Abs[a]^t]^(1/t),
  68.       t===Infinity,   Max[Abs[a]],
  69.       True, Return["ERROR: Norm type not recognized"]
  70.     ],
  71.     MatrixQ[a], Which[
  72.       t===1,          Max[Apply[Plus,Abs[a]]],
  73.       t===2,          Max[SingularValues[N[a]][[2]]],
  74.       t===Infinity,   Max[Apply[Plus,Abs[Transpose[a]]]],
  75.       t===F,          Apply[Plus,Flatten[a]^2]^(1/2),
  76.       True, Return["ERROR: Norm type not recognized"]
  77.     ],
  78.     True, Return["ERROR: Selection must be a vector or matrix"]
  79.   ]
  80.  
  81. tcicondnum[m_] :=
  82.   If[ !(MatrixQ[m] && Dimensions[m][[1]]==Dimensions[m][[2]]), 
  83.     Return["ERROR: selection must be a square matrix"],
  84.     tcinorm[m]*tcinorm[Inverse[m]]
  85.   ]
  86.  
  87. tciSVD[m_] :=
  88.   Module[{mm, uu, dd},
  89.     mm = SingularValues[N[m]];
  90.     uu = Transpose[mm[[1]]];
  91.     dd = DiagonalMatrix[mm[[2]]];
  92.     {uu,dd,mm[[3]]}
  93.   ]
  94.  
  95. tcisvals[m_] := SingularValues[N[m]][[2]];
  96.  
  97. tciQR[m_] :=
  98.   Module[{mm,qq},
  99.     mm = QRDecomposition[N[m]];
  100.     qq = Conjugate[Transpose[mm[[1]]]];
  101.     {qq,mm[[2]]}
  102.   ]
  103.  
  104. lufactor[m_] := LinearAlgebra`GaussianElimination`LUFactor[m];
  105. tciLUfactor[m_] :=
  106.   Module[{mm, ll, pp, uu, nr, nc, xx, yy},
  107.     If[!MatrixQ[m], Return["ERROR: Input not matrix"] ];
  108.     nr = Dimensions[m][[1]];  nc = Dimensions[m][[2]];
  109.     If[nr!=nc, Return["ERROR: Input not square matrix"] ];
  110.     mm = lufactor[m];
  111.     xx = mm[[1]];
  112.     yy = mm[[2]];
  113.     uu = Table[If[j<i,0,xx[[yy[[i]],j]] ],{i,nr},{j,nc}];
  114.     ll = Table[Which[j>i,0,j==i,1,j<i,xx[[yy[[i]],j]] ],{i,nr},{j,nc}];
  115.     pp = Table[If[i==yy[[j]],1,0],{i,nr},{j,nc}];
  116.     {pp, ll, uu}
  117.   ]
  118.  
  119. tciBand[x_,m_,n_] := 
  120.   With[{ k = Length[x], p = Ceiling[Length[x]/2]},
  121.     ff[i_,j_] := Which[j-i+p<1,0, j-i+p>k,0, True,x[[j-i+p]] ];
  122.     Array[ff,{m,n}] ]
  123.  
  124. (* GENERAL *)
  125.  
  126. tcimpy2[a_,b_] := 
  127.   If[ (MatrixQ[a]||VectorQ[a])&&(MatrixQ[b]||VectorQ[b]), a.b, a*b]
  128.  
  129. tcimpy[m_:ListQ] :=
  130.   Module[{ t = m[[1]]},
  131.     Do[ t = tcimpy2[t,m[[i]]], {i,2,Length[m]} ] ;
  132.     t
  133.   ]
  134.  
  135. squarematrix[a_] := MatrixQ[a] && Dimensions[a][[1]]==Dimensions[a][[2]]
  136.  
  137. tciadd[a_,b_] :=
  138.   Module[{c},
  139.     If[ squarematrix[a] && NumberQ[b],
  140.       c = IdentityMatrix[Dimensions[a][[1]] ];
  141.       Return[a+b*c] ];
  142.     If[ squarematrix[b] && NumberQ[a],
  143.       c = IdentityMatrix[Dimensions[b][[1]] ];
  144.       Return[b+a*c] ];
  145.     a+b
  146.   ] 
  147.  
  148. tcipwr[x_,n_,m_:0] :=
  149.   Module[{ w={}, v={}, r=x, nn=n},
  150.     If[ MatrixQ[n], 
  151.       If[ Dimensions[n][[1]] == 1, v = n[[1]]];
  152.       If[ Dimensions[n][[2]] == 1, v = Transpose[n][[1]]]
  153.     ];
  154.     If[ MatrixQ[x] && IntegerQ[n],
  155.       If[ Dimensions[x][[1]] == Dimensions[x][[2]],
  156.         If[m===0,
  157.           Return[MatrixPower[x,n]],
  158.           If[ !PrimeQ[m], Return["ERROR: modulus must be prime"] ];
  159.           If[n < 0, nn=-n; r=Inverse[x, Modulus->m]];
  160.           Return[MatrixPower[r,nn]]                                           
  161.         ],
  162.         Return[x^n]
  163.       ]
  164.     ];
  165.     If[ MatrixQ[x] && (n===H || n===T || n===C),
  166.       If[m===0,
  167.         If[n===T, Return[Transpose[x]] ];
  168.         If[n===C, Return[ComplexExpand[Conjugate[x]]] ];
  169.         If[n===H, Return[ComplexExpand[Conjugate[Transpose[x]]]]],
  170.         Return["ERROR: cannot have modulus"];
  171.       ]
  172.     ];
  173.     If[ MatrixQ[x], 
  174.       If[ Dimensions[x][[1]] == 1, w = x[[1]]];
  175.       If[ Dimensions[x][[2]] == 1, w = Transpose[x][[1]]]
  176.     ];
  177.     If[ (Length[v] > 1) && (Length[w] > 1), 
  178.       Return[ Outer[Power,w,v]],
  179.       Return[x^n]
  180.     ]
  181.   ]
  182.  
  183. tcifactor2[x_,m_:0] := 
  184.   Module[{ yy },
  185.     If[ IntegerQ[x],
  186.       yy = Apply[pow, FactorInteger[x], 2]; Apply[Times,yy],
  187.       If[m===0, Factor[x,Trig->True],Factor[x,Modulus->m]] 
  188.     ]
  189.   ]
  190.  
  191. tcifactor[x_,m_:0] := 
  192.   Which[ MatrixQ[x], Map[tcifactor2[#,m]&,x,2], 
  193.            ListQ[x], Map[tcifactor2[#,m]&,x],
  194.                True, tcifactor2[x,m] 
  195.   ]
  196. (*
  197. tcigetvars[m_] :=
  198.   Module[{ qqv = {}},
  199.     Do[ qqv = Union[qqv, Variables[m[[i,1]]], Variables[m[[i,2]]] ],
  200.               {i,1,Length[m]} ] ;
  201.     Return[qqv]
  202.   ]
  203. *)
  204.  
  205. tcigetvars[x_Symbol] := {x}
  206.  
  207. tcigetvars[x_] := Union @ Cases[Level[x, -1], _Symbol?(Context[#]==="Global`"&)]
  208.  
  209. tcisolve[m_,v_:0] := 
  210.   If[v===0,
  211.     Module[ {qqvars = tcigetvars[m]},
  212.       If[ Length[qqvars] == Length[m],
  213.         Simplify[Solve[m, qqvars]],
  214.         Return["needvars"]
  215.       ]
  216.     ],
  217.     Simplify[Solve[m, v]]
  218.   ]
  219.  
  220. tcigcd2[a_,b_] :=
  221.   If[ IntegerQ[a] && IntegerQ[b], GCD[a,b], PolynomialGCD[a,b] ]
  222.  
  223. tcigcd[a_] := 
  224.   Module[{ t = a[[1]]},
  225.     Do[ t = tcigcd2[t,a[[i]]], {i,2,Length[a]}] ;
  226.     Return[t] ]
  227.  
  228. tcilcm2[a_,b_] :=
  229.   If[ IntegerQ[a] && IntegerQ[b], LCM[a,b], PolynomialLCM[a,b] ]
  230.  
  231. tcilcm[a_] := 
  232.   Module[{ t = a[[1]]},
  233.     Do[ t = tcilcm2[t,a[[i]]], {i,2,Length[a]}] ;
  234.     Return[t] ]
  235.  
  236. tcimixnum[x_] :=
  237.   Module[ {quo, rem, a = Abs[x]},
  238.     quo = Floor[a]; rem = a - quo; quo = Sign[x]*quo;
  239.     If[ quo == 0, Return[rem], Return[mixnum[quo,rem]] ]
  240.   ]
  241.  
  242. tciseries[e_,d_] :=
  243.   Module[ {pow, bb = Series[e,d]},
  244.     pow = bb[[5]]/bb[[6]];
  245.     Return[{Normal[bb] + order[(bb[[1]]-bb[[2]])^pow]} ]
  246.   ]
  247.  
  248. tciODEseries[eqns_, idata_, Y_, x_, x0_, power_] :=
  249.   Module[{ sol, k, vars, U=Y, ee=eqns, id=idata, bb },
  250.     Do[
  251.       U = D[U,x];
  252.       vars = U/.x->x0;
  253.       sol = Solve[(ee/.x->x0)/.id, vars];
  254.       id = Union[ id, sol[[1]] ];
  255.       ee = D[ee, x ],
  256.       {k, 1, power-1}
  257.     ];
  258.     bb = Normal[Series[Y, {x,x0,power-1}]/.id] + order[(x-x0)^power];
  259.     MapThread[#1 == #2 & ,{Y,bb}]
  260.   ]
  261.  
  262. tcifindroot[e_,d_:0] :=
  263.   Module[ {len, ok, var, cc, qq, inf, vlist=tcigetvars[e]},
  264.     If[d===0,
  265.       len = Length[vlist];
  266.       ok = ListQ[e] && Length[e]==len;
  267.       If[!ok, Return["ERROR: system inconsistent"] ];
  268.       If[len==1, FindRoot[e, Append[vlist, 10.*Random[]] ],
  269.         cc = Transpose[{vlist,Table[10.*Random[],{len}]}];
  270.         Apply[FindRoot[e,##]&,cc]
  271.       ],
  272.       ok = ListQ[e] && MatrixQ[d];
  273.       If[!ok, Return["ERROR: data inconsistent"] ];
  274.       ok = Length[e]==Length[d];
  275.       If[!ok, Return["ERROR: constrain all variables or none"] ];
  276.       inf = MemberQ[Flatten[d],Infinity] || MemberQ[Flatten[d],-Infinity];
  277.       If[inf, Return["ERROR: infinite ranges not allowed"] ];
  278.       len = Length[e];
  279.       If[len==1, cc=d; cc[[1,2]]=.5*(cc[[1,3]]+cc[[1,4]]),
  280.         qq = Transpose[d]; qq[[2]]=.5*(qq[[3]]+qq[[4]]);
  281.         cc = Transpose[Take[qq,2]]
  282.       ];
  283.       Apply[FindRoot[e,##]&,cc]
  284.     ]
  285.   ]
  286.  
  287. (* SIMPLEX *)
  288.  
  289. tciminimize[m_,c_] := 
  290.   Module[ { vars },
  291.     vars = Union[Variables[m],tcigetvars[c]];
  292.     ConstrainedMin[m, c, vars]
  293.   ]
  294.  
  295. tcimaximize[m_,c_] := 
  296.   Module[ { vars },
  297.     vars = Union[Variables[m],tcigetvars[c]];
  298.     ConstrainedMax[m, c, vars]
  299.   ]
  300.  
  301. (* POLYNOMIALS *)
  302.  
  303. tcireverse[a_] :=
  304.   Reverse[Apply[List,a]]
  305.  
  306. tcipolydiv[p_,d_,v_:0] :=
  307.   Module[{var},
  308.     If[v===0,
  309.       vv = Variables[p];
  310.       If[Length[vv]==1,var=vv[[1]],Return["needvars"]],
  311.       var=v];
  312.     PolynomialQuotient[p,d,var]+PolynomialRemainder[p,d,var]/d
  313.   ]
  314.  
  315. tciroot[p_,v_:0] :=
  316.   Module[{var,vv,zz},
  317.     If[v===0,
  318.       vv = Variables[p];
  319.       If[Length[vv]==1,var=vv[[1]],Return["needvars"]],
  320.       var=v];
  321.     zz = Roots[p==0,var];
  322.     If[Head[zz]===Roots,
  323.       NSolve[p==0,var],
  324.       Roots[p==0,var] /. lhs_ == rhs_ :> lhs == ComplexExpand[rhs]
  325.     ]
  326.   ]
  327.  
  328. (* CALCULUS *)
  329. tciapart[p_,v_:0] :=
  330.   Module[{var},
  331.     If[v===0,
  332.       vv = Variables[p];
  333.       If[Length[vv]==1,var=vv[[1]],Return["needvars"]],
  334.       var=v];
  335.     Apart[p,var]
  336.   ]
  337.  
  338. tcisum[e_,v_:0,lim_:0] := 
  339.   Module[ {var, f, ans, vlist=Take[Join[{n},tcigetvars[e]],-1]},
  340.     If[v===0, var=vlist[[1]], var=v, ul];
  341.     If[lim===0, lim={1,var}];
  342.     f = e /. var->xx;
  343.     ans = Sum[f, Evaluate[Join[{xx},lim]] ];
  344.     Simplify[ ans /. xx->var ]
  345.   ]
  346.  
  347. tciprod[e_,v_:0,lim_:0] := 
  348.   Module[ {var, f, ans, vlist=Take[Join[{n},tcigetvars[e]],-1]},
  349.     If[v===0, var=vlist[[1]], var=v];
  350.     If[lim===0, lim={1,var}];
  351.     f = e /. var->xx;
  352.     ans = Product[f, Evaluate[Join[{xx},lim]] ];
  353.     Simplify[ ans /. xx->var ]
  354.   ]
  355.  
  356. tciint[e_,v_:0,lim_:0] := 
  357.   Module[ {var, f, ans, vlist=Take[Join[{x},tcigetvars[e]],-1]},
  358.     If[v===0, var=vlist[[1]], var=v];
  359.     If[lim===0, Return[ Integrate[e,var] ] ];
  360.     f = e /. var->xx;
  361.     ans = Integrate[f, Join[{xx},lim] ];
  362.     Simplify[ ans /. xx->var ]
  363.   ]
  364.  
  365. tcintparts[e_,v_,d_] :=
  366.   Module[ {bb = Simplify[e/d], cc, dd, up, low },
  367.     If[ListQ[v], var = v[[1]], var = v];
  368.     cc =Integrate[bb,var];
  369.     dd =D[d,var ];
  370.     If[ListQ[v],
  371.       up = cc*d /.var->v[[3]];
  372.       low = cc*d /.var->v[[2]];
  373.       Return[up - low - int[cc*dd, v]],
  374.       Return[cc*d - int[cc*dd,var ]]
  375.     ]
  376.   ]
  377.  
  378. tciNDSolve[e_,f_,a_,sol_,v_] :=
  379.   Module[ { },
  380.     Clear /@ a;
  381.     sol = NDSolve[e,f,v];
  382.     MapThread[(#1[x_] := Evaluate[#2[x] /. sol])&, {a, f}];
  383.     a
  384.   ]
  385.  
  386.  
  387. tcichangevar[d_,e_,v_] :=
  388.   Module[ {bb , cc, dd, var, newvar=d[[1]], up, low },
  389.     If[ListQ[v], var = v[[1]], var = v];
  390.     dd =d[[2]];
  391.     bb = D[dd,var];
  392.     cc = e/bb /.Solve[d,var];
  393.     If[ListQ[cc], cc = Last[cc]];
  394.     If[ListQ[v],
  395.       up = dd /.var->v[[3]];
  396.       low = dd /.var->v[[2]];
  397.       Return[int[cc, {newvar,low,up}]],
  398.       Return[int[cc, newvar]]
  399.     ]
  400.   ]
  401.  
  402. (* VECTORS *)
  403.  
  404. curl[x_] :=  Calculus`VectorAnalysis`Curl[x]
  405.  
  406. div[x_] :=  Calculus`VectorAnalysis`Div[x]
  407.  
  408. xprod[x_,y_] :=  Calculus`VectorAnalysis`CrossProduct[x,y]
  409.  
  410. tcidprod[a_,b_] := 
  411.   If[ VectorQ[a] && VectorQ[b], a.b, a*b]
  412.  
  413. tcixprod[a_,b_] := 
  414.   If[ VectorQ[a] && VectorQ[b], xprod[a,b], a*b]
  415.  
  416. tcihessian[e_] :=
  417.   Module[ {vars = Variables[e]},
  418.     Array[D[e,vars[[#1]],vars[[#2]]]&,{Length[vars],Length[vars]}]
  419.   ]
  420.  
  421. tcijacobian[a_,v_:0] :=
  422.   Module[{aa, vars, gg},
  423.     If[MatrixQ[a],
  424.       If[Dimensions[a][[2]]==1, aa = Transpose[a][[1]],
  425.         If[Dimensions[a][[1]]==1, aa = a[[1]], 
  426.           Return[ERROR: not row or column matrix] ]
  427.       ], aa = a
  428.     ];
  429.     If[v===0, vars = Variables[aa];
  430.       If[ Length[vars] != Length[aa], Return["needvars"]],
  431.       vars = v
  432.     ];
  433.     gg[i_,j_] := D[aa[[i]],vars[[j]]];
  434.     Array[gg,{Length[vars],Length[vars]}]
  435.   ]
  436.  
  437. tcipotential[a_,v_:0] :=
  438.   Module[{aa, ff, gg, hh, vars, ans, c },
  439.     If[MatrixQ[a],
  440.       If[Dimensions[a][[2]]==1, aa = Transpose[a][[1]],
  441.         If[Dimensions[a][[1]]==1, aa = a[[1]], 
  442.           Return[ERROR: not a row or column matrix] ]
  443.       ], aa = a
  444.     ];
  445.     If[v===0, vars = Variables[aa];
  446.       If[ Length[aa] != 3, Return["ERROR: not a 3D vector field"]];
  447.       If[ Length[vars] != 3, Return["needvars"]],
  448.       vars = v
  449.     ];
  450.     If[tcicurl[aa,vars]!={0,0,0}, Return[0]];
  451.     ans = Integrate[aa[[1]], vars[[1]]];
  452.     Print[ans];
  453.     ans += Integrate[aa[[2]]-D[ans,vars[[2]]], vars[[2]]];
  454.     Print[ans];
  455.     ans += Integrate[aa[[3]]-D[ans,vars[[3]]], vars[[3]]]
  456.   ]
  457.  
  458. tcicurl[a_] := 
  459.   Module[{aa, bb, cc, dd, v },
  460.     If[MatrixQ[a],
  461.       If[Dimensions[a][[2]]==1, aa = Transpose[a][[1]],
  462.         If[Dimensions[a][[1]]==1, aa = a[[1]], 
  463.           Return[ERROR: not a row or column matrix] ]
  464.       ], aa = a
  465.     ];
  466.     curl[aa]
  467.   ]
  468.  
  469. tcidiverge[a_] :=
  470.   Module[{aa},
  471.     If[MatrixQ[a],
  472.       If[Dimensions[a][[2]]==1, aa = Transpose[a][[1]],
  473.         If[Dimensions[a][[1]]==1, aa = a[[1]], 
  474.           Return[ERROR: not a row or column matrix] ]
  475.       ], aa = a
  476.     ];
  477.     div[aa]
  478.   ]
  479.  
  480. (* STATISTICS *)
  481.  
  482. mean[x_] :=  Statistics`DescriptiveStatistics`Mean[x]
  483.  
  484. median[x_] :=  Statistics`DescriptiveStatistics`Median[x]
  485.  
  486. mode[x_] :=  Statistics`DescriptiveStatistics`Mode[x]
  487.  
  488. stddev[x_] :=  Statistics`DescriptiveStatistics`StandardDeviation[x]
  489.  
  490. meandev[x_] :=  Statistics`DescriptiveStatistics`MeanDeviation[x]
  491.  
  492. variance[x_] :=  Statistics`DescriptiveStatistics`Variance[x]
  493.  
  494. tcimean[m_] :=
  495.   If[ MatrixQ[m], 
  496.     If[ NumberQ[ m[[1,1]] ],
  497.       Map[mean,Transpose[m]], 
  498.       Map[mean,Transpose[Drop[m,1]] ]
  499.     ],  
  500.     mean[m]
  501.   ]
  502.  
  503. tcimedian[m_] :=
  504.   If[ MatrixQ[m], 
  505.     If[ NumberQ[ m[[1,1]] ],
  506.       Map[median,Transpose[m]], 
  507.       Map[median,Transpose[Drop[m,1]] ]
  508.     ],  
  509.     median[m]
  510.   ]
  511.  
  512. tcimode[m_] :=
  513.   If[ MatrixQ[m], 
  514.     If[ NumberQ[ m[[1,1]] ],
  515.       Map[mode,Transpose[m]], 
  516.       Map[mode,Transpose[Drop[m,1]] ]
  517.     ],  
  518.     mode[m]
  519.   ]
  520.  
  521. tcistddev[m_] :=
  522.   If[ MatrixQ[m], 
  523.     If[ NumberQ[ m[[1,1]] ],
  524.       Map[stddev,Transpose[m]], 
  525.       Map[stddev,Transpose[Drop[m,1]] ]
  526.     ],  
  527.     stddev[m]
  528.   ]
  529.  
  530.  
  531. tcimeandev[m_] :=
  532.   If[ MatrixQ[m], 
  533.     If[ NumberQ[ m[[1,1]] ],
  534.       Map[meandev,Transpose[m]], 
  535.       Map[meandev,Transpose[Drop[m,1]] ]
  536.     ],  
  537.     meandev[m]
  538.   ]
  539.  
  540. tcivariance[m_] :=
  541.   If[ MatrixQ[m], 
  542.     If[ NumberQ[ m[[1,1]] ],
  543.       Map[variance,Transpose[m]], 
  544.       Map[variance,Transpose[Drop[m,1]] ]
  545.     ],  
  546.     variance[m]
  547.   ]
  548.  
  549. (* Stat functions that don't need the Descriptive Stat package *)
  550.  
  551. tcicorr[m_] := 
  552.   If[ !MatrixQ[m], 
  553.     Return["ERROR: selection must be a matrix"],
  554.     Module[ {nrows, xx, tt, oo, mm = N[m]},
  555.       If[ !NumberQ[ m[[1,1]] ], mm = Drop[m,1], _];
  556.       nrows = Dimensions[mm][[1]];
  557.       tt = Transpose[mm];
  558.       oo = Outer[Times,Table[1,{i,nrows}],Map[mean,tt]];
  559.       xx = (mm - oo).DiagonalMatrix[1/Map[stddev,tt]];
  560.       (Transpose[xx].xx)/(nrows-1)
  561.     ]
  562.   ]
  563.  
  564. tcicovar[m_] := 
  565.   If[ !MatrixQ[m], 
  566.     Return["ERROR: selection must be a matrix"],
  567.     Module[ {nrows, dd, mm = N[m]},
  568.       If[ !NumberQ[ m[[1,1]] ], mm = Drop[m,1]];
  569.       nrows = Dimensions[mm][[1]];
  570.       dd = Map[mean,Transpose[mm]];
  571.       (Transpose[mm].mm)/nrows - Outer[Times,dd,dd]
  572.     ]
  573.   ]
  574.  
  575. tcimoment[m_,p_,a_] :=
  576.   Module[ {nr, nc, bb, cc, mm}, 
  577.     If[MatrixQ[m], mm=m;
  578.       If[ !NumberQ[ m[[1,1]] ], mm = Drop[m,1]];
  579.       nr = Dimensions[mm][[1]]; nc = Dimensions[mm][[2]];
  580.         If[a===M, bb=Map[mean,Transpose[mm]], bb=Table[a,{nc}]];
  581.         cc = (mm-Outer[Times,Table[1,{nr}],bb])^p;
  582.         Return[(Transpose[cc] . Table[1,{nr}])/nr],
  583.       nr = Length[m];
  584.         If[a===M, bb=Table[mean[m],{nr}], bb=Table[a,{nr}]];
  585.         Apply[Plus, ((m-bb)^p)]/nr
  586.     ]
  587.   ]
  588.     
  589. tciquantile[m_,q_] :=
  590.   If[ ((q>1) || (q<0)), 
  591.     Return["ERROR: quantile must lie on [0,1]."],
  592.     If[ MatrixQ[m], 
  593.       If[ NumberQ[ m[[1,1]] ],
  594.         Return[{q,Map[tciquan2[#,q]&,Transpose[m]]}], 
  595.         Return[{q,Map[tciquan2[#,q]&,Transpose[Drop[m,1]]]} ]
  596.       ],  
  597.       {q,tciquan2[m,q]}
  598.     ]
  599.   ]
  600.  
  601. tciquan2[m_,q_] :=
  602.   If[ !VectorQ[m,NumberQ], Return["ERROR: Bad data."],
  603.     Module[ {n = Length[m], ii, ss, mm = Sort[m] },
  604.       ii = Floor[n*q+.5];
  605.       ss = n*q+.5-ii;
  606.       If[ii<1, Return[ mm[[1]] ]];
  607.       If[ii>n, Return[ mm[[n]] ]];
  608.       Return[ mm[[ii]]*(1-ss) + mm[[ii+1]]*ss ]
  609.     ]
  610.   ]
  611.  
  612. tcipolyfit[a_,n_,c_] := 
  613.   If[ !(MatrixQ[a] && Dimensions[a][[2]]==2), 
  614.     Return["ERROR: selection must be a 2 column matrix"],
  615.     Module[{ var, var1, aa = a },
  616.       If[ c===0, aa=Transpose[RotateRight[Transpose[a],1]] ];
  617.       var = aa[[1,1]]; var1 = aa[[1,2]];
  618.       If[ NumberQ[ var ], 
  619.         Return[{y,"=",Fit[a,Table[x^i,{i,0,n}],x] }],
  620.         Return[{var1, "=", Fit[Drop[aa,1],Table[ var^i,{i,0,n}],var] }]
  621.       ]  
  622.     ]
  623.   ]
  624.  
  625. tcimregress[a_,n_,c_] := 
  626.   If[ !MatrixQ[a], Return["ERROR: selection must be a matrix"],
  627.     If[ NumberQ[ a[[1,1]] ], Return["ERROR: columns must be labelled"],
  628.       Module[ {vars, var1, key, aa = a},
  629.         If[c===0, aa=Transpose[RotateLeft[Transpose[a],1]] ];
  630.         vars = Drop[aa[[1]],-1]; var1 = Take[aa[[1]],-1];
  631.         If[n===0, key=vars, key=Join[vars,{1}] ];
  632.         Return[ {var1,"=", Fit[Drop[aa,1],key,vars]} ]
  633.       ]  
  634.     ]
  635.   ]
  636.  
  637. (* PLOTS *)
  638.       
  639. tciview[r_,theta_,phi_] :=
  640.   Module[ {x, y, z},
  641.     x = N[r*Cos[theta*Pi/180.]*Sin[phi*Pi/180.]];
  642.     y = N[r*Sin[theta*Pi/180.]*Sin[phi*Pi/180.]];
  643.     z = N[r*Cos[phi*Pi/180.]];
  644.     {x,y,z}
  645.   ]
  646.  
  647. tciviewlimits[plotlist_]:=
  648.   Module[{outlist={},inlist={},newlist,m,p,j,k,n,approx},
  649.     n=Length[plotlist];
  650.     newlist=Flatten[FullOptions[plotlist,PlotRange]];
  651.     If[Length[newlist]/n==4,m=4,m=6];
  652.     Do[outlist=Append[outlist,Min[Part[newlist,Table[m*(k-1)+p,{k,n}]]]];
  653.        outlist=Append[outlist,Max[Part[newlist,Table[m*(k-1)+p+1,{k,n}]]]],
  654.        {p,1,m-1,2}];
  655.     If[m==4,outlist=Join[outlist,{0,0}]];
  656.     outlist
  657.   ]