I have recently been experimenting with using Axiom for maximum likelihood estimation,
borrowing code (with permission) from "Compact Numerical Methods" by Nash.
\begin{spad}
)abbrev package OPTIM OptimNash
++ Author: M. Clements
++ Date Created: 24 August 2008
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classification:
++ Keywords:
++ References:
++ Description:
++ Numerical optimisation using methods from Nash
++ Permission granted by Nash to translate from his Pascal code to Axiom
axisResults ==> Record(X: Vector Float, Fmin:Float, Lowerfun: Boolean)
nmminResults ==> Record(X:Vector Float, Fmin:Float, Fail:Boolean)
OptimNash: Exports == Implementation where
Exports == with
Rosen: Vector Float -> Float
++ Rosen(Bvec) returns the Rosen function
nmmin: (Vector Float, (Vector Float -> Float), Float) -> nmminResults
++ nmmin(Bvec, fminfn, intol) returns values for the optimised function
axissrch :(Vector Float, Vector Float -> Float) -> axisResults
++ axissrch(Bvec, fminfn) returns the axis search method
demo: () -> nmminResults
++ demo() returns the rsults of a simple demo using the Rosen function
demo2: () -> nmminResults
++ demo2() returns the results of a more complex demo
--optimFun: (Expression Float,List Equation Expression Float) -> nmminResults
Implementation == add
--writeln(Classical starting point (-1.2,1)
)
Rosen(Bvec: Vector Float):Float ==
(Bvec(2)-Bvec(1)^2)^2*100.0+(1.0-Bvec(1))^2
--nmmin(vector([11.0, 1.0]), Rosen$OPTIM, -1.0)$OPTIM
nmmin(Bvec: Vector Float,
fminfn: Vector Float -> Float,
intol: Float):nmminResults ==
Pcol ==> 27
Prow ==> 26
alpha ==> 1.0
beta ==> 0.5
gamma ==> 2.0
Calceps ==> 1.0e-35
big ==> 1.0e35 -- as per R/src/main/optim.c
n : Integer := #Bvec
action : String
C : Integer
calcvert : Boolean
convtol : Float
f : Float
funcount : Integer
H : Integer
i : Integer
j : Integer
L : Integer
notcomp : Boolean
n1 : Integer
oldsize : Float
P : Matrix Float := new(Prow,Pcol,0.0)--[1..Prow,1..Pcol]
shrinkfail: Boolean
size : Float
step : Float
temp : Float
trystep : Float
VH,VL,VN : Float
VR : Float
x : Vector Float := Bvec -- ugly declaration
--writeln(Nash Algorithm 19 version 2 1988-03-17
)
--writeln( Nelder Mead polytope direct search function minimiser
)
--writeln(confile,Nash Algorithm 19 version 2 1988-03-17
)
--writeln(confile, Nelder Mead polytope direct search function minimiser
)
fail := false
(f,notcomp) := (fminfn(Bvec),false)
if notcomp then
--writeln(**** Function cannot be evaluated at initial parameters ****
)
--writeln(confile,
-- **** Function cannot be evaluated at initial parameters ****
)
fail := true
else
--writeln(Function value for initial parameters =
,f)
--writeln(confile,Function value for initial parameters =
,f)
if intol<0.0 then intol := Calceps
funcount := 1
convtol := intol*(abs(f)+intol)
--writeln( Scaled convergence tolerance is
,convtol)
--writeln(confile, Scaled convergence tolerance is
,convtol)
n1 := n+1
C := n+2
P(n1,1) := f
for i in 1..n repeat P(i,1) := Bvec(i)
L := 1
size := 0.0
step := 0.0
for i in 1..n repeat if 0.1abs(Bvec(i))>step then step := 0.1abs(Bvec(i))
--writeln(Stepsize computed as
,step)
--writeln(confile,Stepsize computed as
,step)
for j in 2..n1 repeat
action := "BUILD "
for i in 1..n repeat P(i,j) := Bvec(i)
trystep := step
while P(j-1,j)=Bvec(j-1) repeat
P(j-1,j) := Bvec(j-1)+trystep
trystep := trystep*10.0
size := size+trystep
oldsize := size
calcvert := true
shrinkfail := false
repeat
if calcvert then
for j in 1..n1 repeat
if j~=L then
for i in 1..n repeat Bvec(i) := P(i,j)
(f,notcomp) := (fminfn(Bvec),false)
if notcomp or f>=big then f := big
funcount := funcount+1
P(n1,j) := f
calcvert := false
VL := P(n1,L)
VH := VL
H := L
for j in 1..n1 repeat
if j~=L then
f := P(n1,j)
if fVH then
H := j
VH := f
if VH>VL+convtol then
--str(funcount:5,tstr)
--writeln(action,tstr,' ,VH,
,VL)
--writeln(confile,action,tstr,
,VH,
',VL)
VN := betaVL+(1.0-beta)VH
for i in 1..n repeat
temp := -P(i,H)
for j in 1..n1 repeat temp := temp+P(i,j)
P(i,C) := temp/n
for i in 1..n repeat
Bvec(i) := (1.0+alpha)P(i,C)-alphaP(i,H)
(f,notcomp) := (fminfn(Bvec),false)
if notcomp or f>=big then f := big
funcount := funcount+1
action := "REFLECTION "
VR := f
if VRBvec(i)+(1-gamma)P(i,C)
P(i,C) := Bvec(i)
Bvec(i) := f
(f,notcomp) := (fminfn(Bvec),false)
if notcomp or f>=big then f := big
funcount := funcount+1
if fP(i,H)+betaP(i,C)
(f,notcomp) := (fminfn(Bvec),false)
if notcomp or f>=big then f := big
funcount := funcount+1
if f=VH then
action := "SHRINK "
calcvert := true
size := 0.0
for j in 1..n1 repeat
if j~=L then
for i in 1..n repeat
P(i,j) := beta*(P(i,j)-P(i,L))+P(i,L)
size := size+abs(P(i,j)-P(i,L))
if sizePolytope size measure not decreased in shrink)
--writeln(confile,
-- Polytope size measure not decreased in shrink
)
shrinkfail := true
if ((VH<=VL+convtol) or shrinkfail ) then break -- repeat loop
--writeln(Exiting from Alg19.pas Nelder Mead polytope minimiser
)
--writeln(' ,funcount,
function evaluations used')
--writeln(confile,Exiting from Alg19.pas Nelder Mead polytope minimiser
)
--writeln(confile,' ,funcount,
function evaluations used')
fmin := P(n1,L)
for i in 1..n repeat x(i) := P(i,L)
if shrinkfail then fail := true
return [x, fmin, fail]
axissrch(Bvec: Vector Float, fminfn: Vector Float -> Float): axisResults ==
cradius, eps, f, fplus, step, temp, tilt, fmin: Float
i : Integer
--notcomp : Boolean
--writeln(alg20.pas -- axial search
)
--writeln(confile,alg20.pas -- axial search
)
Calceps : Float := 1.0e-35
big : Float := 1.0e35 -- as per R/src/main/optim.c
n :Integer := #(Bvec)::NonNegativeInteger::Integer
fmin := fminfn(Bvec)
eps := Calceps
eps := sqrt(eps)
--writeln(' Axis':6,' Stepsize ':14,'function + ':14,
-- 'function - :14,
rad. of curv.:14,
tilt')
--writeln(confile,' Axis':6,' Stepsize ':14,'function + ':14,
-- 'function - :14,
rad. of curv.:14,
tilt')
lowerfn := false
for i in 1..n repeat
if (not lowerfn) then
temp := Bvec(i)
step := eps*(abs(temp)+eps)
Bvec(i) := temp+step
f := fminfn(Bvec)
if (f>=big) then f := big
--write(i:5,' ,step:12,
,f:12,
)
--write(confile,i:5,
,step:12,
,f:12,
)
if f<fmin then lowerfn := true
if (not lowerfn) then
fplus := f
Bvec(i) := temp-step
f := fminfn(Bvec)
if f>=big then f := big
--write(f:12,
) write(confile,f:12,
')
if f temp := 0.5(fplus-f)/step
fplus := 0.5(fplus+f-2.0fmin)/(stepstep)
if fplus~=0.0 then
cradius := 1.0+temptemp
cradius := cradiussqrt(cradius)/fplus
else
cradius := big
tilt := 45.0*atan(temp)/atan(1.0)
--write(cradius:12,' ,tilt:12)
--write(confile,cradius:12,
',tilt:12)
--writeln writeln(confile)
[Bvec, f, lowerfn]
demo(): nmminResults ==
nmmin([-1.2, 1.0], Rosen, -1.0)
demo2(): nmminResults ==
--banner:=dr1920.pas -- driver for Nelder-Mead minimisation
--startup
--fminset(n,B,Workdata) {sets up problem and defines starting
-- values of B}
lowerfn:Boolean := false --{safety setting}
B:Vector Float := [-1.2,1.0] --initial values
mynmminResults :nmminResults
myaxisResults : axisResults
--n :Integer := #(B)::NonNegativeInteger::Integer
repeat
mytol:=-1.0 --Note: set the tolerance negative to indicate that
--procedure must obtain an appropriate value.
mynmminResults := nmmin(B,Rosen,mytol) --{minimise the function}
--writeln
--writeln(confile)
--writeln( Minimum function value found =
,Fmin)
--writeln( At parameters
)
--writeln(confile, Minimum function value found =
,Fmin)
--writeln(confile, At parameters
)
--for i in 1..n repeat
--{
--writeln(' B[',i,]=
,X[i])
--writeln(confile,' B[',i,]=
,X[i])
--} --{loop to write out parameters}
B := mynmminResults.X
myaxisResults := axissrch(B, Rosen) --{alg20.pas}
lowerfn := myaxisResults.Lowerfun
if lowerfn then
B := myaxisResults.X
--writeln(Lower function value found
)
--writeln(confile,Lower function value found
)
if (not lowerfn) then break
return mynmminResults
--flush(confile) close(confile)
--if infname<>con
then close(infile)
--{dr1920.pas -- Nelder Mead minimisation with axial search}
-- optimFun(e:Expression Float,initial:List Equation Expression Float):nmminResults ==
-- vars := [lhs(x) for x in initial]
-- vals := [rhs(x) for x in initial]::Vector Float
-- nmmin(vals, (x +-> eval(e, vars, x)), -1)
\end{spad}
\begin{spad}
)abbrev package SOMESTAT SomeStatisticsPackage
SomeStatisticsPackage: Exports == Implementation where
DF ==> DoubleFloat
Exports ==> with
pgamma:(Float,Float) -> Float
++ pgamma(y,p) returns the incomplete gamma function (Alg AS 147)
erf:(Float) -> Float
++ erf(x) returns the error function
erfc:(Float) -> Float
++ erfc(x) returns the complementary error function
pnorm:(Float) -> Float
++ pnorm(x) returns the CDF for the normal(0,1) function
pnorm:(Float,Float,Float) -> Float
++ pnorm(x,mean,sd) returns the CDF for a normal(mean,sd) function
pchisq:(Float,Float) -> Float
ppois:(Integer,Float) -> Float
ppois2:(Integer,Float) -> Float
rpois:(Float) -> Integer
++ rpois(lambda) returns a random Poisson variable with mean lambda
rbinom:(Integer,Float) -> Integer
++ rbinom(n,p) returns a random binomial variable from n trials
++ and probability p
mean:(Vector Float) -> Float
mean:(List Float) -> Float
mean:(Vector DF) -> DF
mean:(List DF) -> DF
var:(Vector Float) -> Float
Beta:(Float,Float) -> Float
betai:(Float,Float,Float) -> Float
betacf:(Float,Float,Float) -> Float
pbeta:(Float,Float,Float) -> Float
++ pbeta(x,a,b) returns the incomplete beta function betai(a,b,x)
pt:(Float,Float) -> Float
++ pt(x,df) returns the CDF for a t distribution
pbinom:(Float,Float,Float) -> Float
++ pbinom(k,n,p) returns the CDF for 0..k events given n trials with probability p
pf:(Float,Float,Float) -> Float
++ pf(f,df1,df2) returns the CDF of the F distribution with df1 and df1 degrees of freedom
Implementation ==> add
import FloatSpecialFunctions
import RandomFloatDistributions
pgamma(y,p) ==
eps : Float := (10.0^(-digits()$Float+1))$Float
g : Float := 0.0
if (y=0.0 or p=0.0) then return g
if (y < 0.0 or p < 0.0) then error "Invalid arguments"
a : Float := p + 1.0
f : Float := exp(p log(y) - logGamma(a) - y)
if (f < eps) then
return g
c : Float := 1.0
g : Float := 1.0
a : Float := p
while (c > eps g) repeat
a := a+1.0
c := c(y / a)
g := g+c
g := gf
return g
erf(x:Float):Float == if x<0.0 then -erf(-x) else pgamma(x^2,0.5)
erfc(x:Float):Float == 1-erf(x)
pnorm(x:Float):Float == 0.5(1+erf(x/sqrt(2.0)))
pnorm(x:Float,mu:Float,sigma:Float):Float == 0.5(1+erf((x-mu)/sigma/sqrt(2.0)))
pchisq(x:Float,df:Float):Float == pgamma(x/2,df/2)
ppois(y:Integer,lambda:Float):Float ==
if y<0 then return 0.0
reduce(_+,[exp(-lambda)lambda^yi/factorial(yi) for yi in 0..y])$List(Float)
ppois2(y:Integer,lambda:Float):Float == 1-pgamma(lambda,y::Float+1)
rpois(lambda:Float):Integer ==
cumP : Float := 0
x : Float := uniform01()
for i in 0.. repeat
cumP := cumP + exp(-lambda)lambda^i/factorial(i)
if x1.0 then
error "bad argument x in betai"
bt := if x=0.0 or x=1.0 then 0.0 else
x^a(1.0-x)^b/Beta(a,b)
-- exp(logGamma(a+b)$FloatSpecialFunctions-
-- logGamma(a)$FloatSpecialFunctions-
-- logGamma(b)$FloatSpecialFunctions+
-- alog(x)+blog(1.0-x))
if x<(a+1.0)/(a+b+2.0) then return bt betacf(a,b,x)/a else
return 1.0-btbetacf(b,a,1.0-x)/b
betacf(a:Float,b:Float,x:Float):Float ==
itmax : Integer := 1000
eps : Float := (10.0^(-digits()$Float+1))$Float
qab : Float := a+b
qap : Float := a+1.0
qam : Float := a-1.0
c : Float := 1.0
d : Float := 1.0-qabx/qap
d := 1.0/d
h : Float := d
for m in 1..itmax repeat
em : Float := m::Float
m2 := 2.0em
aa : Float := m(b-em)x/((qam+m2)(a+m2))
d := 1.0+aad
c := 1.0+aa/c
d := 1.0/d
h := hdc
aa := -(a+em)(qab+em)x/((a+m2)(qap+m2))
d := 1.0+aad
c := 1.0+aa/c
d := 1.0/d
del : Float := dc
h := hdel
if abs(del-1.0)f))
\end{spad}
Now, let's try this.
\begin{axiom}
nmminResults ==> Record(X:Vector Float, Fmin:Float, Fail:Boolean)
optimFun(e:Expression Float,initial:List Equation Expression Float):nmminResults ==
vars := [lhs(x) for x in initial]
vals := [rhs(x) for x in initial]::Vector Float
nmmin(vals, ((x:Vector Float):Float +-> eval(e, vars, x)), -1)$OPTIM
-- some examples
dpois := exp(-lambda)*lambda^y/Gamma(y+1)
negll1 := -log(eval(dpois, y=10.0))
optimFun(negll1, [lambda=5.0])
optimFun((b-a^2)^2*100.0+(1.0-a)^2, [a=11.0,b=1.0])
x := [xi::Float for xi in 1..5]
Y := [10, 15, 20, 25, 31]
negll1 := reduce(+, [eval(-log(dpois),[y=Y(i),lambda=exp(alpha+beta*x(i))]) for i in 1..#Y]);
out1 := optimFun(negll1, [alpha=2, beta=0])
-- R: summary(glm(c(10, 15, 20, 25, 31) ~ I(1:5), family=poisson))
mleResult ==> Record(coef:Matrix Float,vars:List Symbol, AIC:Float)
mle(negll:Expression Float, initial:List Equation Expression Float):mleResult ==
vars := [lhs(x) for x in initial]::List Symbol
vals := [rhs(x) for x in initial]::Vector Float
out : nmminResults := nmmin(vals, (x +-> eval(negll, vars, x)), -1)$OPTIM
estimate := out.X
covMatrix := inverse(matrix([[eval(D(D(negll,x),y),vars,estimate) for y in vars] for x in vars]))
se := [sqrt(covMatrix(i,i)) for i in 1..nrows(covMatrix)]
zscore := [estimate(i)/se(i) for i in 1..nrows(covMatrix)]
pvalues := [2(1-pnorm(zscore(i)::Float)$SOMESTAT) for i in 1..nrows(covMatrix)]
AIC := out.Fmin2+2*#(out.X)
coef := transpose(matrix([estimate,se,zscore,pvalues]))
return [coef,vars,AIC]
mle(negll1, [alpha=2.0, beta=0.0])
negll2 := -log(eval(dpois, y=10.0))
mle(negll2, [lambda=5.0])
\end{axiom}
As an aside: can we calculate an elasticity from logistic regression?
\begin{axiom}
Y:=exp(alpha+betakPk)/(1+exp(alpha+betakPk))
D(Y,Pk)/Y
D(Y,Pk)/Y/(1-Y)
\end{axiom}
Some or all expressions may not have rendered properly,
because Axiom returned the following error:
Error: export FRICAS=/usr/local/lib/fricas/target/x86_64-unknown-linux; export ALDORROOT=/usr/local/aldor/linux/1.1.0; export PATH=$ALDORROOT/bin:$PATH; export HOME=/var/zope2/var/LatexWiki; ulimit -t 600; export LD_LIBRARY_PATH=/usr/local/lib/fricas/target/x86_64-unknown-linux/lib; LANG=en_US.UTF-8 $FRICAS/bin/FRICASsys < /var/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/7245897146845666469-25px.axm
/bin/sh: /usr/local/lib/fricas/target/x86_64-unknown-linux/bin/FRICASsys: not found
Some or all expressions may not have rendered properly,
because Latex returned the following error:
This is pdfTeXk, Version 3.141592-1.40.3 (Web2C 7.5.6)
\write18 enabled.
%&-line parsing enabled.
entering extended mode
(./3789055917543086017-16.0px.tex
LaTeX2e <2005/12/01>
Babel <v3.8h> and hyphenation patterns for english, usenglishmax, dumylang, noh
yphenation, arabic, farsi, croatian, ukrainian, russian, bulgarian, czech, slov
ak, danish, dutch, finnish, basque, french, german, ngerman, ibycus, greek, mon
ogreek, ancientgreek, hungarian, italian, latin, mongolian, norsk, icelandic, i
nterlingua, turkish, coptic, romanian, welsh, serbian, slovenian, estonian, esp
eranto, uppersorbian, indonesian, polish, portuguese, spanish, catalan, galicia
n, swedish, ukenglish, pinyin, loaded.
(/usr/share/texmf-texlive/tex/latex/base/article.cls
Document Class: article 2005/09/16 v1.4f Standard LaTeX document class
(/usr/share/texmf-texlive/tex/latex/base/size12.clo))
(/usr/share/texmf-texlive/tex/latex/ucs/ucs.sty
(/usr/share/texmf-texlive/tex/latex/ucs/data/uni-global.def))
(/usr/share/texmf-texlive/tex/latex/base/inputenc.sty
(/usr/share/texmf-texlive/tex/latex/ucs/utf8x.def))
(/usr/share/texmf-texlive/tex/latex/bbm/bbm.sty)
(/usr/share/texmf-texlive/tex/latex/jknapltx/mathrsfs.sty)
(/usr/share/texmf-texlive/tex/latex/base/fontenc.sty
(/usr/share/texmf-texlive/tex/latex/base/t1enc.def))
(/usr/share/texmf-texlive/tex/latex/pstricks/pstricks.sty
(/usr/share/texmf-texlive/tex/generic/pstricks/pstricks.tex
`PSTricks' v1.15 <2006/12/22> (tvz)
(/usr/share/texmf-texlive/tex/generic/pstricks/pstricks.con))
(/usr/share/texmf/tex/latex/xcolor/xcolor.sty
(/etc/texmf/tex/latex/config/color.cfg)
(/usr/share/texmf-texlive/tex/latex/graphics/dvips.def)
(/usr/share/texmf-texlive/tex/latex/graphics/dvipsnam.def)))
(/usr/share/texmf-texlive/tex/latex/graphics/epsfig.sty
(/usr/share/texmf-texlive/tex/latex/graphics/graphicx.sty
(/usr/share/texmf-texlive/tex/latex/graphics/keyval.sty)
(/usr/share/texmf-texlive/tex/latex/graphics/graphics.sty
(/usr/share/texmf-texlive/tex/latex/graphics/trig.sty)
(/etc/texmf/tex/latex/config/graphics.cfg))))
(/usr/share/texmf-texlive/tex/latex/pst-grad/pst-grad.sty
(/usr/share/texmf-texlive/tex/generic/pst-grad/pst-grad.tex
(/usr/share/texmf-texlive/tex/latex/xkeyval/pst-xkey.tex
(/usr/share/texmf-texlive/tex/latex/xkeyval/xkeyval.sty
(/usr/share/texmf-texlive/tex/latex/xkeyval/xkeyval.tex)))
`pst-plot' v1.05, 2006/11/04 (tvz,dg,hv)))
(/usr/share/texmf-texlive/tex/latex/pstricks/pst-plot.sty
(/usr/share/texmf-texlive/tex/generic/pstricks/pst-plot.tex
v97 patch 2, 1999/12/12
(/usr/share/texmf-texlive/tex/generic/multido/multido.tex
v1.41, 2004/05/18 <tvz>)))
(/usr/share/texmf-texlive/tex/latex/geometry/geometry.sty
(/usr/share/texmf-texlive/tex/xelatex/xetexconfig/geometry.cfg)
Package geometry Warning: `lmargin' and `rmargin' result in NEGATIVE (-108.405p
t).
`width' should be shortened in length.
) (/usr/share/texmf-texlive/tex/latex/amsmath/amsmath.sty
For additional information on amsmath, use the `? option.
(/usr/share/texmf-texlive/tex/latex/amsmath/amstext.sty
(/usr/share/texmf-texlive/tex/latex/amsmath/amsgen.sty))
(/usr/share/texmf-texlive/tex/latex/amsmath/amsbsy.sty)
(/usr/share/texmf-texlive/tex/latex/amsmath/amsopn.sty))
(/usr/share/texmf-texlive/tex/latex/amsfonts/amsfonts.sty)
(/usr/share/texmf-texlive/tex/latex/amsfonts/amssymb.sty)
(/usr/share/texmf-texlive/tex/latex/amscls/amsthm.sty)
(/usr/share/texmf-texlive/tex/latex/setspace/setspace.sty
Package: `setspace
6.7 <2000/12/01>
) (/usr/share/texmf-texlive/tex/generic/xypic/xy.sty
(/usr/share/texmf-texlive/tex/generic/xypic/xy.tex Bootstrap'ing: catcodes,
docmode, (/usr/share/texmf-texlive/tex/generic/xypic/xyrecat.tex)
(/usr/share/texmf-texlive/tex/generic/xypic/xyidioms.tex)
Xy-pic version 3.7 <1999/02/16>
Copyright (c) 1991-1998 by Kristoffer H. Rose <krisrose@ens-lyon.fr>
Xy-pic is free software: see the User's Guide for details.
Loading kernel: messages; fonts; allocations: state, direction,
utility macros; pictures: \xy, positions, objects, decorations;
kernel objects: directionals, circles, text; options; algorithms: directions,
edges, connections; Xy-pic loaded))
(/usr/share/texmf-texlive/tex/generic/xypic/xyall.tex
Xy-pic option: All features v.3.3
(/usr/share/texmf-texlive/tex/generic/xypic/xycurve.tex
Xy-pic option: Curve and Spline extension v.3.7 curve, circles, loaded)
(/usr/share/texmf-texlive/tex/generic/xypic/xyframe.tex
Xy-pic option: Frame and Bracket extension v.3.7 loaded)
(/usr/share/texmf-texlive/tex/generic/xypic/xycmtip.tex
Xy-pic option: Computer Modern tip extension v.3.3
(/usr/share/texmf-texlive/tex/generic/xypic/xytips.tex
Xy-pic option: More Tips extension v.3.3 loaded) loaded)
(/usr/share/texmf-texlive/tex/generic/xypic/xyline.tex
Xy-pic option: Line styles extension v.3.6 loaded)
(/usr/share/texmf-texlive/tex/generic/xypic/xyrotate.tex
Xy-pic option: Rotate and Scale extension v.3.3 loaded)
(/usr/share/texmf-texlive/tex/generic/xypic/xycolor.tex
Xy-pic option: Colour extension v.3.3 loaded)
(/usr/share/texmf-texlive/tex/generic/xypic/xymatrix.tex
Xy-pic option: Matrix feature v.3.4 loaded)
(/usr/share/texmf-texlive/tex/generic/xypic/xyarrow.tex
Xy-pic option: Arrow and Path feature v.3.5 path, \ar, loaded)
(/usr/share/texmf-texlive/tex/generic/xypic/xygraph.tex
Xy-pic option: Graph feature v.3.7 loaded) loaded)
(/usr/share/texmf-texlive/tex/latex/tools/verbatim.sty)
(/usr/share/texmf/tex/latex/graphviz/graphviz.sty
(/usr/share/texmf-texlive/tex/latex/psfrag/psfrag.sty))
(/usr/share/texmf/tex/latex/sagetex.sty
Writing sage input file 3789055917543086017-16.0px.sage
) (/usr/share/texmf-texlive/tex/latex/gnuplottex/gnuplottex.sty
(/usr/share/texmf-texlive/tex/latex/base/latexsym.sty)
(/usr/share/texmf-texlive/tex/latex/moreverb/moreverb.sty)
(/usr/share/texmf-texlive/tex/latex/base/ifthen.sty))
(./3789055917543086017-16.0px.aux)
(/usr/share/texmf-texlive/tex/latex/ucs/ucsencs.def)
Missing $ inserted.
<inserted text>
$
l.153 (Bvec(2)-Bvec(1)^
2)^2*100.0+(1.0-Bvec(1))^2
(/usr/share/texmf-texlive/tex/latex/jknapltx/ursfs.fd)
(/usr/share/texmf-texlive/tex/latex/amsfonts/umsa.fd)
(/usr/share/texmf-texlive/tex/latex/amsfonts/umsb.fd)
(/usr/share/texmf-texlive/tex/latex/base/ulasy.fd)
Missing $ inserted.
<inserted text>
$
l.154
You can't use `macro parameter character #' in horizontal mode.
l.166 n : Integer := #
Bvec
Overfull \hbox (27.39969pt too wide) in paragraph at lines 155--223
[]\T1/cmr/m/n/12 --nmmin(vector([11.0, 1.0]), Rosen$\OML/cmm/m/it/12 OPTIM; \OM
S/cmsy/m/n/12 ^^@\OT1/cmr/m/n/12 1\OML/cmm/m/it/12 :\OT1/cmr/m/n/12 0)$\T1/cmr/
m/n/12 OPTIM nm-min(Bvec: Vec-tor Float, fminfn: Vec-tor Float -> Float, in-tol
: Float):nmminResults
[1]
You can't use `macro parameter character #' in horizontal mode.
l.333 n :Integer := #
(Bvec)::NonNegativeInteger::Integer
You can't use `macro parameter character #' in horizontal mode.
l.388 --n :Integer := #
(B)::NonNegativeInteger::Integer
[2]
Missing $ inserted.
<inserted text>
$
l.464 eps : Float := (10.0^
(-digits()$Float+1))$Float
Missing $ inserted.
<inserted text>
$
l.494 cumP := cumP + exp(-lambda)*lambda^
i/factorial(i)
Missing $ inserted.
<inserted text>
$
l.504 mean(x:Vector Float):Float == reduce(_
+,x)/#x
You can't use `macro parameter character #' in math mode.
l.504 ...n(x:Vector Float):Float == reduce(_+,x)/#
x
You can't use `macro parameter character #' in math mode.
l.505 mean(x:List Float):Float == reduce(_+,x)/#
x
You can't use `macro parameter character #' in math mode.
l.506 mean(x:Vector DF):DF == reduce(_+,x)/#
x
You can't use `macro parameter character #' in math mode.
l.507 mean(x:List DF):DF == reduce(_+,x)/#
x
You can't use `macro parameter character #' in math mode.
l.510 reduce(_+,[(x(i)-meanx)^2 for i in 1..#
x])$List(Float)/(#x-1)
You can't use `macro parameter character #' in horizontal mode.
l.510 ...-meanx)^2 for i in 1..#x])$List(Float)/(#
x-1)
Missing $ inserted.
<inserted text>
$
l.527 eps : Float := (10.0^
(-digits()$Float+1))$Float
Missing $ inserted.
<inserted text>
$
l.557 \end{spad}
\newpage
Overfull \hbox (11.6583pt too wide) in paragraph at lines 422--557
\OML/cmm/m/it/12 Float \OT1/cmr/m/n/12 == \OML/cmm/m/it/12 ify < \OT1/cmr/m/n/1
2 0\OML/cmm/m/it/12 thenreturn\OT1/cmr/m/n/12 0\OML/cmm/m/it/12 :\OT1/cmr/m/n/1
2 0\OML/cmm/m/it/12 reduce\OT1/cmr/m/n/12 ([]\OML/cmm/m/it/12 ; \OT1/cmr/m/n/12
[\OML/cmm/m/it/12 exp\OT1/cmr/m/n/12 (\OMS/cmsy/m/n/12 ^^@\OML/cmm/m/it/12 lam
bda\OT1/cmr/m/n/12 ) \OMS/cmsy/m/n/12 ^^C \OML/cmm/m/it/12 lambda[]i=factorial\
OT1/cmr/m/n/12 (\OML/cmm/m/it/12 yi\OT1/cmr/m/n/12 )\OML/cmm/m/it/12 foryiin\OT
1/cmr/m/n/12 0\OML/cmm/m/it/12 ::y\OT1/cmr/m/n/12 ])$\T1/cmr/m/n/12 List(Float)
ppois2(y:Integer,lambda:Float):Float
Overfull \hbox (65.58339pt too wide) in paragraph at lines 422--557
\OML/cmm/m/it/12 meanx \OT1/cmr/m/n/12 : \OML/cmm/m/it/12 Float \OT1/cmr/m/n/12
:= \OML/cmm/m/it/12 mean\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 x\OT1/cmr/m/n/12 )\O
ML/cmm/m/it/12 reduce\OT1/cmr/m/n/12 ([]\OML/cmm/m/it/12 ; \OT1/cmr/m/n/12 [(\O
ML/cmm/m/it/12 x\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 i\OT1/cmr/m/n/12 ) \OMS/cmsy/
m/n/12 ^^@ \OML/cmm/m/it/12 meanx\OT1/cmr/m/n/12 )[]\OML/cmm/m/it/12 foriin\OT1
/cmr/m/n/12 1\OML/cmm/m/it/12 ::x\OT1/cmr/m/n/12 ])$\T1/cmr/m/n/12 List(Float)/
(x-1) --Gamma(x:Float):Float == Gamma(x)$\OML/cmm/m/it/12 FloatSpecialFunctions
Beta\OT1/cmr/m/n/12 (\OML/cmm/m/it/12 a \OT1/cmr/m/n/12 :
[3] (/usr/share/texmf-texlive/tex/latex/base/t1cmtt.fd)
LaTeX Warning: Characters dropped after `\end{axiom}' on input line 593.
LaTeX Warning: Characters dropped after `\end{axiom}' on input line 598.
[4] [5] (./3789055917543086017-16.0px.aux) )
(see the transcript file for additional information)
Output written on 3789055917543086017-16.0px.dvi (5 pages, 21476 bytes).
Transcript written on 3789055917543086017-16.0px.log.