|
|
|
last edited 14 years ago by Mark Clements |
Edit detail for SandBoxMLE revision 2 of 3
| 1 2 3 | ||
|
Editor: Mark Clements
Time: 2008/12/04 18:42:35 GMT-8 |
||
| Note: | ||
added: \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 := g*f 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 x<cumP then return i rbinom(size:Integer,p:Float):Integer == y : Integer := 0 for i in 1..size repeat if uniform01()<p then y := y + 1 return(y) -- rBernoulli := [(if uniform01()<p then 1::Integer else 0::Integer) -- for i in 1..size] -- reduce(_+,rBernoulli)$List(Integer) mean(x:Vector Float):Float == reduce(_+,x)/#x mean(x:List Float):Float == reduce(_+,x)/#x mean(x:Vector DF):DF == reduce(_+,x)/#x mean(x:List DF):DF == reduce(_+,x)/#x var(x:Vector Float):Float == meanx : Float := mean(x) reduce(_+,[(x(i)-meanx)**2 for i in 1..#x])$List(Float)/(#x-1) --Gamma(x:Float):Float == Gamma(x)$FloatSpecialFunctions Beta(a:Float,b:Float):Float == Gamma(a)*Gamma(b)/Gamma(a+b) -- see Press et al (1992) betai(a:Float,b:Float,x:Float):Float == if x<0.0 or x>1.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+ -- a*log(x)+b*log(1.0-x)) if x<(a+1.0)/(a+b+2.0) then return bt * betacf(a,b,x)/a else return 1.0-bt*betacf(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-qab*x/qap d := 1.0/d h : Float := d for m in 1..itmax repeat em : Float := m::Float m2 := 2.0*em aa : Float := m*(b-em)*x/((qam+m2)*(a+m2)) d := 1.0+aa*d c := 1.0+aa/c d := 1.0/d h := h*d*c aa := -(a+em)*(qab+em)*x/((a+m2)*(qap+m2)) d := 1.0+aa*d c := 1.0+aa/c d := 1.0/d del : Float := d*c h := h*del if abs(del-1.0)<eps then return h error "WARNING: a or b too big, or itmax too small" --return h pbeta(x:Float,a:Float,b:Float):Float == betai(a,b,x) pt(x:Float,df:Float):Float == 1-betai(df/2.0,0.5,df/(df+x**2))/2.0 pbinom(k:Float,n:Float,p:Float):Float == 1-betai(k+1.0,n-k,p) pf(f:Float,df1:Float,df2:Float):Float == 1-betai(df2/2.0,df1/2.0,df2/(df2+df1*f)) \end{spad} removed: -VF ==> Vector Float changed: -optimFun(e:Expression Float,initial:List Equation Expression Float):nmminResults == optimFun(e:Expression Float,initial:List Equation Expression Float):nmminResults == added: -- some examples dpois := exp(-lambda)*lambda^y/Gamma(y+1) negll1 := -log(eval(dpois, y=10.0)) optimFun(negll1, [lambda=5.0]) added: 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.Fmin*2+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])
I have recently been experimenting with using Axiom for maximum likelihood estimation, borrowing code (with permission) from "Compact Numerical Methods" by Nash.
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.1*abs(Bvec(i))>step then step := 0.1*abs(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 f<VL then
L := j
VL := f
if f>VH 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 := beta*VL+(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)-alpha*P(i,H)
(f,notcomp) := (fminfn(Bvec),false)
if notcomp or f>=big then f := big
funcount := funcount+1
action := "REFLECTION "
VR := f
if VR<VL then
P(n1,C) := f
for i in 1..n repeat
f := gamma*Bvec(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 f<VR then
for i in 1..n repeat P(i,H) := Bvec(i)
P(n1,H) := f
action := "EXTENSION "
else
for i in 1..n repeat P(i,H) := P(i,C)
P(n1,H) := VR
else
action := "HI-REDUCTION "
if VR<VH then
for i in 1..n repeat P(i,H) := Bvec(i)
P(n1,H) := VR
action := "LO-REDUCTION "
for i in 1..n repeat Bvec(i) := (1-beta)*P(i,H)+beta*P(i,C)
(f,notcomp) := (fminfn(Bvec),false)
if notcomp or f>=big then f := big
funcount := funcount+1
if f<P(n1,H) then
for i in 1..n repeat P(i,H) := Bvec(i)
P(n1,H) := f
else
if VR>=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 size<oldsize then
shrinkfail := false
oldsize := size
else
--writeln('Polytope 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<fmin then lowerfn := true
if (not lowerfn) then
Bvec(i) := temp
temp := 0.5*(fplus-f)/step
fplus := 0.5*(fplus+f-2.0*fmin)/(step*step)
if fplus~=0.0 then
cradius := 1.0+temp*temp
cradius := cradius*sqrt(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)
spad
Compiling FriCAS source code from file
/var/zope2/var/LatexWiki/3925042998800764367-25px001.spad using
old system compiler.
OPTIM abbreviates package OptimNash
processing macro definition axisResults ==> Record(X: Vector Float,Fmin: Float,Lowerfun: Boolean)
processing macro definition nmminResults ==> Record(X: Vector Float,Fmin: Float,Fail: Boolean)
------------------------------------------------------------------------
initializing NRLIB OPTIM for OptimNash
compiling into NRLIB OPTIM
compiling exported Rosen : Vector Float -> Float
****** comp fails at level 5 with expression: ******
error in function Rosen
(+
(* (** (- (|Bvec| 2) | << | (** (|Bvec| 1) 2) | >> |) 2)
((|elt| (|Float|) |float|) 100 0 10))
(** (- ((|elt| (|Float|) |float|) 1 0 10) (|Bvec| 1)) 2))
****** level 5 ******
$x:= (** (Bvec (One)) 2)
$m:= $EmptyMode
$f:=
((((|Bvec| # #) (|#| #) (* #) (+ #) ...)))
>> Apparent user error:
Cannot coerce 2
of mode (PositiveInteger)
to mode (UniversalSegment (Integer))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 := g*f
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 x<cumP then return i
rbinom(size:Integer,p:Float):Integer ==
y : Integer := 0
for i in 1..size repeat
if uniform01()<p then y := y + 1
return(y)
-- rBernoulli := [(if uniform01()<p then 1::Integer else 0::Integer)
-- for i in 1..size]
-- reduce(_+,rBernoulli)$List(Integer)
mean(x:Vector Float):Float == reduce(_+,x)/#x
mean(x:List Float):Float == reduce(_+,x)/#x
mean(x:Vector DF):DF == reduce(_+,x)/#x
mean(x:List DF):DF == reduce(_+,x)/#x
var(x:Vector Float):Float ==
meanx : Float := mean(x)
reduce(_+,[(x(i)-meanx)**2 for i in 1..#x])$List(Float)/(#x-1)
--Gamma(x:Float):Float == Gamma(x)$FloatSpecialFunctions
Beta(a:Float,b:Float):Float == Gamma(a)*Gamma(b)/Gamma(a+b)
-- see Press et al (1992)
betai(a:Float,b:Float,x:Float):Float ==
if x<0.0 or x>1.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+
-- a*log(x)+b*log(1.0-x))
if x<(a+1.0)/(a+b+2.0) then return bt * betacf(a,b,x)/a else
return 1.0-bt*betacf(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-qab*x/qap
d := 1.0/d
h : Float := d
for m in 1..itmax repeat
em : Float := m::Float
m2 := 2.0*em
aa : Float := m*(b-em)*x/((qam+m2)*(a+m2))
d := 1.0+aa*d
c := 1.0+aa/c
d := 1.0/d
h := h*d*c
aa := -(a+em)*(qab+em)*x/((a+m2)*(qap+m2))
d := 1.0+aa*d
c := 1.0+aa/c
d := 1.0/d
del : Float := d*c
h := h*del
if abs(del-1.0)<eps then return h
error "WARNING: a or b too big, or itmax too small"
--return h
pbeta(x:Float,a:Float,b:Float):Float == betai(a,b,x)
pt(x:Float,df:Float):Float == 1-betai(df/2.0,0.5,df/(df+x**2))/2.0
pbinom(k:Float,n:Float,p:Float):Float == 1-betai(k+1.0,n-k,p)
pf(f:Float,df1:Float,df2:Float):Float ==
1-betai(df2/2.0,df1/2.0,df2/(df2+df1*f))
spad
Compiling FriCAS source code from file
/var/zope2/var/LatexWiki/4835601963657463339-25px002.spad using
old system compiler.
SOMESTAT abbreviates package SomeStatisticsPackage
processing macro definition DF ==> DoubleFloat
processing macro definition Exports ==> -- the constructor category
processing macro definition Implementation ==> -- the constructor capsule
------------------------------------------------------------------------
initializing NRLIB SOMESTAT for SomeStatisticsPackage
compiling into NRLIB SOMESTAT
importing FloatSpecialFunctions
importing RandomFloatDistributions
compiling exported pgamma : (Float,Float) -> Float
****** comp fails at level 3 with expression: ******
error in function pgamma
(SEQ
(LET (|:| |eps| (|Float|))
| << |
((|elt| (|Float|) **) ((|elt| (|Float|) |float|) 10 0 10)
(+ (- ((|elt| (|Float|) |digits|))) 1))
| >> |)
(LET (|:| |g| (|Float|))
((|elt| (|Float|) |float|) 0 0 10))
(SEQ
(LET #1=#:G725
(= |y| ((|elt| (|Float|) |float|) 0 0 10)))
(|exit| 1
(IF #1#
(|return| 1 |g|)
(SEQ
(LET #2=#:G726
(= |p| ((|elt| (|Float|) |float|) 0 0 10)))
(|exit| 1
(IF #2#
(|return| 1 |g|)
|noBranch|))))))
(SEQ
(LET #3=#:G727
(< |y| ((|elt| (|Float|) |float|) 0 0 10)))
(|exit| 1
(IF #3#
(|error| #4="Invalid arguments")
(SEQ
(LET #5=#:G728
(< |p| ((|elt| (|Float|) |float|) 0 0 10)))
(|exit| 1
(IF #5#
(|error| #4#)
|noBranch|))))))
(LET (|:| |a| (|Float|))
(+ |p| ((|elt| (|Float|) |float|) 1 0 10)))
(LET (|:| |f| (|Float|))
(|exp| (- (- (* |p| (|log| |y|)) (|logGamma| |a|)) |y|)))
(IF (< |f| |eps|)
(|return| 1 |g|)
|noBranch|)
(LET (|:| |c| (|Float|))
((|elt| (|Float|) |float|) 1 0 10))
(LET (|:| |g| (|Float|))
((|elt| (|Float|) |float|) 1 0 10))
(LET (|:| |a| (|Float|))
|p|)
(REPEAT (WHILE (< (* |eps| |g|) |c|))
(SEQ
(LET |a|
(+ |a| ((|elt| (|Float|) |float|) 1 0 10)))
(LET |c|
(* |c| (/ |y| |a|)))
(|exit| 1
(LET |g|
(+ |g| |c|)))))
(LET |g|
(* |g| |f|))
(|exit| 1 (|return| 1 |g|)))
****** level 3 ******
$x:= ((elt (Float) **) ((elt (Float) float) 10 (Zero) 10) (+ (- ((elt (Float) digits))) (One)))
$m:= $EmptyMode
$f:=
((((|eps| #) (|p| # . #1=#) (|y| # #) (|p| . #1#) ...)))
>> Apparent user error:
NoValueMode
is an unknown modeNow, let's try this.
axiom
nmminResults ==> Record(X:Vector Float, Fmin:Float, Fail:Boolean)
Type: Void
axiom
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)$OPTIM
Function declaration optimFun : (Expression(Float),List(Equation(
Expression(Float)))) -> Record(X: Vector(Float),Fmin: Float,Fail
: Boolean) has been added to workspace.
Type: Void
axiom
-- some examples dpois := exp(-lambda)*lambda^y/Gamma(y+1)
 | (1) |
Type: Expression(Integer)
axiom
negll1 := -log(eval(dpois, y=10.0))
 | (2) |
Type: Expression(Float)
axiom
optimFun(negll1, [lambda=5.0])
You cannot now use OptimNash in the context you have it. optimFun((b-a**2)**2*100.0+(1.0-a)**2, [a=11.0,b=1.0])
>> System error: The function |*2;optimFun;1;initial| is undefined.
As an aside: can we calculate an elasticity from logistic regression?
axiom
Y:=exp(alpha+betak*Pk)/(1+exp(alpha+betak*Pk))
 | (3) |
Type: Expression(Integer)
axiom
D(Y,Pk)/Y
 | (4) |
Type: Expression(Integer)
axiom
D(Y,Pk)/Y/(1-Y)
 | (5) |
Type: Expression(Integer)