[q-lang-cvs] qcalc/examples eds-example.qcalc,NONE,1.1

[Albert G. <ag...@us...>]

Update of /cvsroot/q-lang/qcalc/examples
In directory sc8-pr-cvs16.sourceforge.net:/tmp/cvs-serv30963

Added Files:
	eds-example.qcalc 
Log Message:
added Eddie Rucker's examples

--- NEW FILE: eds-example.qcalc ---
// qcalc 1.0, created Mon Nov 26 22:17:41 2007 -*-Q-*- -*- coding: UTF-8 -*-
// [((0,0),"Take the gamma of B1:"),((0,1),"4.3"),((0,2),"gamma B1 ="),((0,3),"= gamma 4.3"),((1,0),"Factorial is based on gamma:"),((1,1),"8"),((1,2),"factorial B2 ="),((1,3),"= factorial B2"),((3,0),"How many ways can we choose"),((3,1),"Size of group:"),((3,2),"8"),((3,3),"Answer:"),((3,4),"= permutations C4 C5"),((4,0),"some elements from a set when order"),((4,1),"Number to choose:"),((4,2),"4"),((5,0),"is important."),((7,0),"How many ways can we choose"),((7,1),"Size of group:"),((7,2),"8"),((7,3),"Answer:"),((7,4),"= combinations C8 C9"),((8,0),"some elements from a set when order"),((8,1),"Number to choose:"),((8,2),"4"),((9,0),"is *not* important."),((11,0),"Binomial Distribution"),((11,1),"Times coin tossed:"),((11,2),"3"),((11,3),"Answer:"),((11,4),"= binomialDist C12 C13 C14"),((12,0),"A coin is tossed C12 times. What is"),((12,1),"Number of Heads:"),((12,2),"2"),((13,0),"the probability of getting exactly "),((13,1),"Prob. of a Head:"),((13,2),".5"),((14,0),"C13 heads"),((16,0),"Dates"),((16,1),"Values"),((16,2),"Percent Change"),((17,0),"2000"),((17,1),"10"),((18,0),"2001"),((18,1),"15"),((18,2),"= 100*percentChange B18 B19"),((19,0),"2002"),((19,1),"12"),((19,2),"= 100*percentChange B19 B20"),((20,0),"2003"),((20,1),"16"),((20,2),"= 100*percentChange B20 B21"),((21,0),"2004"),((21,1),"19"),((21,2),"= 100*percentChange B21 B22"),((23,0),"Given a set of numbers"),((23,1),"Numbers"),((24,0),"we can find the following:"),((24,1),"3"),((25,1),"2"),((26,1),"4"),((27,1),"5"),((28,1),"3"),((29,1),"4"),((30,1),"2"),((31,1),"3"),((32,1),"2"),((33,0),"Mean:"),((33,1),"= mean B25:B33"),((34,0),"Mode:"),((34,1),"= mode (sort (<) B25:B33)"),((35,0),"Median:"),((35,1),"= median B25:B33"),((36,0),"Lower Quartile:"),((36,1),"= lowerQuartile B25:B33"),((37,0),"Upper Quartile:"),((37,1),"= upperQuartile B25:B33"),((38,0),"Range:"),((38,1),"= range B25:B33"),((39,0),"Population Variance:"),((39,1),"= populationVariance B25:B33"),((40,0),"Population Standard Deviation:"),((40,1),"= populationSTD B25:B33"),((41,0),"Sample Variance:"),((41,1),"= sampleVariance B25:B33"),((42,0),"Sample Standard Deviation:"),((42,1),"= sampleSTD B25:B33"),((44,0),"Find outliers of a group"),((44,1),"Numbers"),((45,1),"1"),((46,1),"30"),((47,1),"3"),((48,1),"-20"),((49,1),"5"),((50,1),"6"),((51,1),"7"),((52,0),"outliers = "),((52,1),"= outliers (sort (<) B46:B52)"),((54,0),"Regression"),((54,1),"X"),((54,2),"Y"),((55,0),"Find the best fitting line for the"),((55,1),"1"),((55,2),"4"),((56,0),"points to the right."),((56,1),"2"),((56,2),"6"),((57,1),"3"),((57,2),"7"),((58,1),"4"),((58,2),"9"),((59,1),"5"),((59,2),"12"),((60,1),"7"),((60,2),"13"),((62,0),"Least Squares:"),((62,1),"= regression B56:B61 C56:C61"),((63,0),"Value of regression at point B63"),((63,1),"9"),((63,2),"= regression B56:B61 C56:C61 B64"),((64,0),"correlation of coefficient"),((64,1),"= correlation B56:B61 C56:C61"),((65,0),"SSE"),((65,1),"= _SSE B56:B61 C56:C61"),((66,0),"MSE"),((66,1),"= _MSE B56:B61 C56:C61"),((68,0),"Assuming a Normal distribution"),((69,0),"Mean:"),((69,1),"30"),((70,0),"Standard Deviation:"),((70,1),"5"),((71,1),"23"),((71,2),"41"),((72,0),"Probability X is between B71 and C71"),((72,1),"= normal (zScore B72 B70 B71) (zScore C72 B70 B71)"),((73,0),"Probability X is above B71"),((73,1),"= normalAbove (zScore B72 B70 C72)"),((74,0),"Probability X is below B71"),((74,1),"= normalBelow (zScore B72 B70 B71)")]
// []
// [(0,380),(1,181),(2,164),(3,117),(4,117)]
// Start of script. Please do not remove this line.
/* Written by Eddie Rucker 2006-2007
   This library contains various functions that I have
   used for various projects and some I thought up for fun.
   I don't promise they will be efficient or be useful and 
   I don't provide a warranty of any kind. 
   
   Use at your own risk!
*/

public var const pi = 3.1415926535897931;

public (\\) X Y @ (-);
public medianSorted, median, lowerQuartileSorted, lowerQuartile, maximum,
  minimum, upperGamma, lowerGamma, percentChange, outliers, deleteOutliers,
  mean, mode, range, sampleVariance, populationVariance, sampleSTD, 
  populationSTD,  correlation, integrate, regression, _SXY, _SSE, _MSE, zScore,
  normal, normalBelow, normalAbove, gamma, factorial, permutations, combinations, 
  binomialDist, beta, erf, efc;

/* remove elements of Y from X (like set difference)
   example [1,2,3,4,5,1,6] \\ [1,3] => [2,4,5,6]
*/

[] \\ Y:List		= [];
Xs:List \\ []		= Xs;
Xs:List \\ [Y|Ys]	= (del_ Xs Y []) \\ Ys;

/* remove an element K from a list */
private del_;
del_ [] K Acc		= reverse Acc;
del_ [K|Ls] K Acc 	= del_ Ls K Acc;
del_ [L|Ls] K Acc 	= del_ Ls K [L|Acc];

/* gamma */
gamma 1			= 1;
gamma Z:Int		= _FAIL_			  if Z <= 0;
  			= (Z-1) * gamma (Z-1)		  otherwise;
gamma Z:Num		= pi / (sin (pi*Z) * gamma (1-Z)) if re Z < 0.5;
	= 2.506628274631 * T^(Z-0.5) * exp(-T)*(0.99999999999980993 + sum [P!N/(Z+N) : N in [0..7]])
		where 
			T = Z + 6.5, 
			P = [676.5203681218851, -1259.1392167224028, 
			     771.32342877765313, -176.61502916214059, 
			     12.507343278686905, -0.13857109526572012, 
			     9.9843695780195716e-6, 1.5056327351493116e-7];

/* factorial of integer N */
factorial N:Int	= gamma (N+1);

/* permutation R from N elements */
permutations N:Int R:Int	
	= _FAIL_ 		 if (N < 0) or else (N < R);
	= permutation_ N (N-R) 1 otherwise;
permutation_ N N S	= S;
permutation_ N R S	= permutation_ (N-1) R (N*S);

/* combination R from N elements */
combinations N:Int R:Int	
	= _FAIL_					if (N < 0) or else (N < R);
	= 0						if N < R;
	= (permutations N Diff) div (factorial Diff) 	if R >= Diff where Diff = N-R;
	= (permutations N R) div (factorial R)		otherwise;

/* binomial distribution with probability P */
binomialDist N:Int X:Int P:Float
	= _FAIL_ 					if (X < 0) or else ((P <= 0.0) and (P >= 1.0));
	= (combinations N X) * (P^X) * (1.0-P)^(N-X)	otherwise;

/* percentage increase/decrease from X to Y */
percentChange X Y = (Y - X)/X;

/* median assuming that L is sorted
   an average of the two medians is returned in case #L is even
*/
medianSorted L:List 
	= _FAIL_		if #L < 3;
where M = (#L) div 2:
	= (L!(M-1) + L!M) / 2	if #L mod 2 = 0;
	= L!M			otherwise;

/* return the median assuming that L is not sorted */
median 	= medianSorted . sort (<);

/* return the mode assuming that L is sorted */
mode [X|Xs] 
where L = reverse (mode_ Xs X 0 0 []):
	= [] if L = [X|Xs];
	= L;

mode_ [] CurrElem OldCount CurrCount Acc
	= Acc            if CurrCount < OldCount;
	= [CurrElem|Acc] if CurrCount = OldCount;
	= [CurrElem]     if CurrCount > OldCount;
mode_ [X|Xs] CurrElem OldCount CurrCount Acc
	= mode_ Xs X OldCount (CurrCount + 1) Acc if X = CurrElem;
	= mode_ Xs X CurrCount 0 [CurrElem]       if CurrCount > OldCount;
	= mode_ Xs X CurrCount 0 [CurrElem|Acc]   if CurrCount = OldCount;
	= mode_ Xs X OldCount 0 Acc               if CurrCount < OldCount;

/* lower quartile assuming that L is sorted */
lowerQuartileSorted L:List = medianSorted (take (#L div 2) L);

/* lower quartile assuming that L is not sorted */
lowerQuartile = lowerQuartileSorted . sort (<);

/* upper quartile assuming that L is sorted */
upperQuartileSorted L:List = medianSorted (drop (ceil (#L+1) div 2) L);

/* upper quartile assuming that L is not sorted */
upperQuartile = upperQuartileSorted . sort (<);

/* outliers assuming that L is sorted
   uses the 1.5*Interquartile method
*/

outliers L:List = filter (\Y. (Y < LB) or else (Y > UB)) L
	where 
		(UQ, LQ) = (upperQuartileSorted L, lowerQuartileSorted L),
		IRQ = 1.5 * (UQ - LQ),
		(LB, UB) = (LQ - IRQ, UQ + IRQ);

/* removes outliers from L; assume that L is sorted */
deleteOutliers L:List = L \\ (outliers L);

/* returns mean of a list L */
mean L:List = (sum L) / #L;

/* maximum of a list */
minimum L:List = foldl1 min L;

/* minimum of a list */
maximum L:List = foldl1 max L;

/* returns the range of list L */
range L:List = (maximum L) - (minimum L);

private sqr;
sqr X = X * X;

/* returns the sample variance of L */
sampleVariance L:List
	= ((sum (map sqr L)) - (sqr (sum L))/N) / (N - 1) where N = #L;

/* sample standard deviation */
sampleSTD L:List = sqrt (sampleVariance L);

/* returns the population variance of L */
populationVariance L:List
	= ((sum (map sqr L)) - (sqr (sum L))/N) / N where N = #L;

/* population standard of L */
populationSTD L:List = sqrt (populationVariance L);

/* returns the covariance of lists X and Y */
covariance X:List Y:List
	= ((sum (zipwith (*) X Y)) - (sum X) * (sum Y) / #X) / (#X - 1);

/* returns the correlation coefficient of Lists X and Y */
correlation X:List Y:List
	= (covariance X Y) / (sqrt (sampleVariance X) * sqrt (sampleVariance Y));

/* Numeric integration using the Recursive Simpson's Rule */

integrate F A:Real B:Real 
	= rsr F A FA B FB 0.0000005 Est
		where 
			FA = F A, 
			FB = F B,
			M = 0.5*(A+B),
			Est = ((abs (M-A)/6.0) * (FA + 4.0*(F M) + FB));

/* rsr = Recursive Simpson's Rule */
private rsr;
rsr F A FA B FB Err Sum
	where 
		M = 0.5*(A+B),
		C = abs (M-A)/6.0,
		FM = F M,
		L = C * (FA + 4.0*F (0.5*(A+M)) + FM),
		R = C * (FM + 4.0*F (0.5*(M+B)) + FB),
		C2 = L+R-Sum:
	= L + R + C2/15.0 if abs C2 <= 15.0*Err;
	= (rsr F A FA M FM (0.5*Err) L) + (rsr F M FM B FB (0.5*Err) R);

/* examples of integration 

   These methods are not the best algorithms! 
   DON'T USE THEM FOR SERIOUS WORK!

   example: lowerGamma 1 10 => 0.999954669331989
   octave: gammai(1,10) => ans = 0.99995

*/

lowerGamma A X 
	= 1.0 - exp(-X) if A = 1;
	= 1.0/gamma A * integrate (\T. exp(-T) * T^E) 0.0 X where E = A - 1;

/* incomplete upper gamma */
upperGamma A X = (gamma A) - lowerGamma A X;

/* error function */
erf X = 1.12837916709551 * integrate (\Z . exp(-Z*Z)) 0.0 X;

/* complementary error function */
erfc X = abs (1.0 - erf X);

/* cumulative chi-square distribution 
   Please don't use this! 

chiSquareCDF X D
  = erf (sqrt (0.5*X)) if D = 1; 
  = (lowerGamma (D/2) (X/2))/gamma(D/2);
*/

/* Beta function */
beta X Y = (gamma X * gamma Y) / gamma (X + Y);


/* regression stuff

   regression X Y => \X1 . A+B*X1

   examples: 
     regression [1,2,3,4] [3,4,5,6] => \X1 . 2.0+1.0*X1
     regression [1,2,3,4] [3,4,5,6] 5 => 7.0
*/

/* sum of squares */
_SXY X:List Y:List = (sum $ zipwith (*) X Y) - (sum X * sum Y) / (#X);

private _SR;
_SR X:List Y:List = (_SXY X Y)^2 / _SXY X X;

/* sum of squares error */
_SSE X:List Y:List = (_SXY Y Y) - _SR X Y;

/* sum of mean squares error */
_MSE X:List Y:List = _SSE X Y / (#X - 2);

regression X:List Y:List
	= \Z . A + B * Z
		where	
			B = (_SXY X Y) / (_SXY X X),
			A = (mean Y) - B * (mean X);


/* Normal Distribution

   normal A B  => probability of P(A<X<B)
   normalBelow A => P(X<A)
   normalAbove A => P(X>A)
*/

zScore X Mean Std = (X-Mean) / Std;

/* normal PDF function */

private normalPDF;
normalPDF Z:Real = 0.3989422804014327 * exp (-0.5*Z*Z);

/* normal A B => probability of P(A<X<B) */ 
normal = integrate normalPDF;

normalBelow Z:Real
	= 0.5 + integrate normalPDF 0.0 Z if Z > 0.0;
	= 0.5 - integrate normalPDF 0.0 Z if Z < 0.0;
	= 0.5;

normalAbove Z:Real = 1.0 - normalBelow Z;