[Axiom-developer] [statistical functions] Old hacks

[unknown]

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
                if S has RadicalCategory then
                        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  
  Implementation == add
        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
                        if S has RadicalCategory 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)
                                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
                        if S has RadicalCategory 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 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
                if S has RadicalCategory then
                        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

  Implementation == add
        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
                        if S has RadicalCategory then
                                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

  Implementation == add
        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
                                
                        if S has RadicalCategory then
                                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}

--
forwarded from http://page.axiom-developer.org/zope/mathaction/address@hidden