axiom-developer
[Top][All Lists]

## [Axiom-developer] [statistical functions] Old hacks

 From: unknown Subject: [Axiom-developer] [statistical functions] Old hacks Date: Sat, 11 Jun 2005 07:02:10 -0500

Changes
http://page.axiom-developer.org/zope/mathaction/StatisticalFunctions/diff
--
Here is some old hack.
\begin{axiom}
)abbrev package STAT StatPackage
++ Description:
++ This package exports statistics utilities
StatPackage(S,A) : Exports == Implementation where
S: SetCategory
A: Collection(S) with (finiteAggregate)
Exports == with
if A has FiniteLinearAggregate(S) and S has OrderedSet then
median: A -> S
++ median(a) median of a collection
if S has Field then
mean: A -> S
++ mean(a) arithmetic mean of a collection
hmean: A -> S
++ hmean(a) harmonic mean of a collection
moment: (A,Integer) -> S
++ moment(a,i) i nth moment of a collection
variance: A -> S
++ variance(a) variance of a collection
center: A -> A
++ center(a) center collection
if A has shallowlyMutable then
center_!: A -> A
++ center!(a) center collection
gmean: A -> S
++ gmean(a) geometric mean of a collection
stdev: A -> S
++ stdev(a) root mean square of a collection
stdize: A -> A
++ stdize(a) standardize a collection
if A has shallowlyMutable then
stdize_!: A -> A
++ stdize!(a) standardize a collection
if A has FiniteLinearAggregate(S) and S has OrderedSet then
median(m)==
n := sort(m)
i:= divide(#m,2)
if (zero? i.remainder) then
j:= elt(n,i.quotient+minIndex(m)-1)
else
j:= elt(n,i.quotient+minIndex(m))
j

if S has Field then
if A has OneDimensionalArrayAggregate(S) then
sums: A ->List S
sums(m)==
j:=0::S
j2:=0::S
for i in minIndex(m)..maxIndex(m) repeat
h:= elt(m,i)
j := j + h
j2 := j2 + h*h
[j,j2]
variance(m)==
i := center(m)
s := sums(i)
k := (s.2 - s.1/#m::S)/(#m::S - 1::S)
k
stdize(m)==
i := mean(m)
j := map(#1-i,m)
s := sums(j)
k := sqrt((s.2 - s.1/#m::S)/(#m::S -
1::S))
map((#1-i)/k,m)
if A has shallowlyMutable then
stdize_!(m)==
i := mean(m)
j := map(#1-i,m)
s := sums(j)
k := sqrt((s.2 -
s.1/#m::S)/(#m::S - 1::S))
map_!((#1-i)/k,m)

else
variance(m)==
i := center(m)
j := reduce("+",i)
j2 := reduce("+", map(#1*#1,i))
k := (j2 - j/#m::S)/(#m::S - 1::S)
k
stdize(m)==
me:= mean(m)
i := map(#1-me,m)
j := reduce("+",i)
j2 := reduce("+", map(#1*#1,i))
k := sqrt((j2 - j/#m::S)/(#m::S - 1::S))
map((#1-me)/k,m)
if A has shallowlyMutable then
stdize_!(m)==
me:= mean(m)
i := map(#1-me,m)
j := reduce("+",i)
j2 := reduce("+", map(#1*#1,i))
k := sqrt((j2 - j/#m::S)/(#m::S
- 1::S))
map_!((#1-me)/k,m)

if A has FiniteLinearAggregate(S) and S has OrderedSet then
median(m)==
n := sort(m)
min := minIndex(m)
i:= divide(#m,2)
if (zero? i.remainder) then
j:=
(elt(n,i.quotient+min-1)+elt(n,i.quotient+min))/2::S
else
j:=elt(n,i.quotient+min)
j
mean(m)==
i := reduce("+", m)
(i / (#m::S))

hmean(m)==
i := map(1::S/#1,m)
j := reduce("+", i)
((#m::S) / j)

moment(m,n)==
n = 0 => 1
n = 1 => mean(m)
i := map(#1**n, m)
j := reduce("+",i)
(j / (#m::S))

center(m)==
i := mean(m)
map(#1-i,m)

if A has shallowlyMutable then
center_!(m)==
i := mean(m)
map_!(#1-i,m)
m
gmean(m)==
i := reduce("*", m)
nthRoot(i,#m)
stdev(m)==
sqrt(variance(m))

)abbrev package IDSTAT IntegralDomainStatPackage
++ Description:
++ This package exports statistics utilities over IntegralDomain
IntegralDomainStatPackage(S,A) : Exports == Implementation where
S: IntegralDomain
A: Collection(S) with (finiteAggregate)
Exports == with
mean: A -> Fraction(S)
++ mean(a) arithmetic mean of a collection
moment: (A,NonNegativeInteger) -> Fraction(S)
++ moment(m,n) nth moment of a collection

mean(m)==
i := reduce("+", m)
i / #m::S
moment(m,n)==
n = 0 => 1
n = 1 => mean(m)
i := map(#1**n, m)
j := reduce("+",i)
j / #m::S
)abbrev package MSTAT MatrixStatPackage
++ Description:
++ This package exports statistics utilities on two dimensional arrays
++ Compute with column major order
MatrixStatPackage(S,Row,Col,M) : Exports == Implementation where
S: SetCategory
Row : FiniteLinearAggregate S
Col : FiniteLinearAggregate S
M   : TwoDimensionalArrayCategory(S,Row,Col)
Exports == with
if Col has shallowlyMutable then
if S has OrderedSet then
median: M -> Col
++ median(a) median of a two dimensional array
if S has Field then
mean: M -> Col
++ mean(a) arithmetic mean of a two dimensional array
hmean: M -> Col
++ hmean(a) harmonic mean of a two dimensional array
variance: M -> Col
++ variance(a) variance of a two dimensional array
moment: (M,Integer) -> Col
++ moment(a,i) i nth moment of a two dimensional array
center: M -> M
++ center(a) center two dimensional array
center_!: M -> M
++ center!(a) center two dimensional array
if M has MatrixCategory(S,Row,Col) then
cov:  M -> M
++ cov(a): covariance matrix
gmean: M -> Col
++ gmean(a) geometric mean of a two dimensional
array
stdev: M -> Col
++ stdev(a) root mean square of a two
dimensional array
stdize: M -> M
++ stdize(a) standardize two dimensional array
stdize_!: M -> M
++ stdize!(a) standardize two dimensional array
if M has MatrixCategory(S,Row,Col) then
cor:  M -> M
++ cor(a): correlation matrix

import Lisp
import StatPackage(S,Col)
if Col has shallowlyMutable then
colReduce: (Col -> S, M) -> Col
colReduce(fx,m)==
ret : Col := new(ncols(m),NIL$Lisp) j := minColIndex(m) - 1 for i in j+1..maxColIndex(m) repeat setelt(ret,i-j,fx(column(m,i))) ret if S has OrderedSet then median(m) == ret := colReduce(median$StatPackage(S,Col),m)
ret
if S has Field then
mean(m)==
ret := colReduce(mean$StatPackage(S,Col),m) ret hmean(m)== ret := colReduce(hmean$StatPackage(S,Col),m)
ret
variance(m)==
ret := colReduce(variance$StatPackage(S,Col),m) ret moment(m,pow)== ret : Col := new(ncols(m),NIL$Lisp)
j := minColIndex(m) - 1
for i in j+1..maxColIndex(m) repeat

setelt(ret,i-j,moment(column(m,i),pow)$StatPackage(S,Col)) ret center(m)== ret : M := new(nrows(m),ncols(m),NIL$Lisp)
j := minColIndex(m) - 1
for i in j+1..maxColIndex(m) repeat

setColumn_!(ret,i-j,center_!(column(m,i))$StatPackage(S,Col)) ret center_!(m)== j := minColIndex(m) - 1 for i in j+1..maxColIndex(m) repeat setColumn_!(m,i-j,center_!(column(m,i))$StatPackage(S,Col))
m
if M has MatrixCategory(S,Row,Col) then
cov(m)==
ret := center(m)
transpose(ret/(nrows(m)-1)::S) * ret

gmean(m)==
colReduce(gmean$StatPackage(S,Col),m) stdev(m)== colReduce(stdev$StatPackage(S,Col),m)
stdize(m)==
ret : M :=
new(nrows(m),ncols(m),NIL$Lisp) j := minColIndex(m) - 1 for i in j+1..maxColIndex(m) repeat setColumn_!(ret,i-j,stdize_!(column(m,i))$StatPackage(S,Col))
ret
stdize_!(m)==
j := minColIndex(m) - 1
for i in j+1..maxColIndex(m) repeat

setColumn_!(m,i-j,stdize_!(column(m,i))\$StatPackage(S,Col))
m
if M has MatrixCategory(S,Row,Col) then
cor(m)==
ret := stdize(m)

(transpose(ret/((nrows(m)-1)::S)) * ret)

\end{axiom}

--