login  home  contents  what's new  discussion  bug reports     help  links  subscribe  changes  refresh  edit

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/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/4784896474666810238-25px001.spad
      using old system compiler.
   OPTIM abbreviates package OptimNash 
******** Spad syntax error detected ********
The prior line was:
37> intol: Float):nmminResults ==
The current line is:
38> Pcol ==> 27
The number of valid tokens is 2. The prior token was #S(TOKEN :SYMBOL == :TYPE KEYWORD :NONBLANK 37) The current token is #S(TOKEN :SYMBOL |;| :TYPE KEYWORD :NONBLANK 37) The next token is #S(TOKEN :SYMBOL |Pcol| :TYPE IDENTIFIER :NONBLANK 38)

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/lib/zope2.10/instance/axiom-wiki/var/LatexWiki/3087042336776025296-25px002.spad
      using old system compiler.
   SOMESTAT abbreviates package SomeStatisticsPackage 
------------------------------------------------------------------------
   initializing NRLIB SOMESTAT for SomeStatisticsPackage 
   compiling into NRLIB SOMESTAT 
   importing FloatSpecialFunctions
   importing RandomFloatDistributions
   compiling exported pgamma : (Float,Float) -> Float
Time: 0 SEC.
compiling exported erf : Float -> Float Time: 0.01 SEC.
compiling exported erfc : Float -> Float Time: 0 SEC.
compiling exported pnorm : Float -> Float Time: 0.02 SEC.
compiling exported pnorm : (Float,Float,Float) -> Float Time: 0 SEC.
compiling exported pchisq : (Float,Float) -> Float Time: 0.01 SEC.
compiling exported ppois : (Integer,Float) -> Float Time: 0.01 SEC.
compiling exported ppois2 : (Integer,Float) -> Float Time: 0.01 SEC.
compiling exported rpois : Float -> Integer Time: 0 SEC.
compiling exported rbinom : (Integer,Float) -> Integer Time: 0 SEC.
compiling exported mean : Vector Float -> Float Time: 0 SEC.
compiling exported mean : List Float -> Float Time: 0 SEC.
compiling exported mean : Vector DoubleFloat -> DoubleFloat Time: 0 SEC.
compiling exported mean : List DoubleFloat -> DoubleFloat Time: 0 SEC.
compiling exported var : Vector Float -> Float Time: 0.01 SEC.
compiling exported Beta : (Float,Float) -> Float Time: 0 SEC.
compiling exported betai : (Float,Float,Float) -> Float Time: 0.02 SEC.
compiling exported betacf : (Float,Float,Float) -> Float Time: 0.02 SEC.
compiling exported pbeta : (Float,Float,Float) -> Float Time: 0.01 SEC.
compiling exported pt : (Float,Float) -> Float Time: 0 SEC.
compiling exported pbinom : (Float,Float,Float) -> Float Time: 0 SEC.
compiling exported pf : (Float,Float,Float) -> Float Time: 0 SEC.
(time taken in buildFunctor: 0)
;;; *** |SomeStatisticsPackage| REDEFINED
;;; *** |SomeStatisticsPackage| REDEFINED Time: 0 SEC.
Warnings: [1] rbinom: y has no value
Cumulative Statistics for Constructor SomeStatisticsPackage Time: 0.12 seconds
finalizing NRLIB SOMESTAT Processing SomeStatisticsPackage for Browser database: --------constructor--------- --------(pgamma ((Float) (Float) (Float)))--------- --------(erf ((Float) (Float)))--------- --------(erfc ((Float) (Float)))--------- --------(pnorm ((Float) (Float)))--------- --------(pnorm ((Float) (Float) (Float) (Float)))--------- --->-->SomeStatisticsPackage((pchisq ((Float) (Float) (Float)))): Not documented!!!! --->-->SomeStatisticsPackage((ppois ((Float) (Integer) (Float)))): Not documented!!!! --->-->SomeStatisticsPackage((ppois2 ((Float) (Integer) (Float)))): Not documented!!!! --------(rpois ((Integer) (Float)))--------- --------(rbinom ((Integer) (Integer) (Float)))--------- --->-->SomeStatisticsPackage((mean ((Float) (Vector (Float))))): Not documented!!!! --->-->SomeStatisticsPackage((mean ((Float) (List (Float))))): Not documented!!!! --->-->SomeStatisticsPackage((mean ((DoubleFloat) (Vector (DoubleFloat))))): Not documented!!!! --->-->SomeStatisticsPackage((mean ((DoubleFloat) (List (DoubleFloat))))): Not documented!!!! --->-->SomeStatisticsPackage((var ((Float) (Vector (Float))))): Not documented!!!! --->-->SomeStatisticsPackage((Beta ((Float) (Float) (Float)))): Not documented!!!! --->-->SomeStatisticsPackage((betai ((Float) (Float) (Float) (Float)))): Not documented!!!! --->-->SomeStatisticsPackage((betacf ((Float) (Float) (Float) (Float)))): Not documented!!!! --------(pbeta ((Float) (Float) (Float) (Float)))--------- --------(pt ((Float) (Float) (Float)))--------- --------(pbinom ((Float) (Float) (Float) (Float)))--------- --------(pf ((Float) (Float) (Float) (Float)))--------- ; compiling file "/var/aw/var/LatexWiki/SOMESTAT.NRLIB/SOMESTAT.lsp" (written 31 JUL 2013 03:50:40 PM):
; /var/aw/var/LatexWiki/SOMESTAT.NRLIB/SOMESTAT.fasl written ; compilation finished in 0:00:00.231 ------------------------------------------------------------------------ SomeStatisticsPackage is now explicitly exposed in frame initial SomeStatisticsPackage will be automatically loaded when needed from /var/aw/var/LatexWiki/SOMESTAT.NRLIB/SOMESTAT

Now, let's try this.

fricas
nmminResults ==> Record(X:Vector Float, Fmin:Float, Fail:Boolean)
Type: Void
fricas
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
Function declaration optimFun : (Expression(Float),List(Equation( Expression(Float)))) -> Record(X: Vector(Float),Fmin: Float,Fail : Boolean) has been added to workspace.
Type: Void
fricas
-- some examples
dpois := exp(-lambda)*lambda^y/Gamma(y+1)

\label{eq1}{{{e}^{- lambda}}\ {{lambda}^{y}}}\over{\Gamma \left({y + 1}\right)}(1)
Type: Expression(Integer)
fricas
negll1 := -log(eval(dpois, y=10.0))

\label{eq2}\begin{array}{@{}l}
\displaystyle
-{{1.0}\ {\log{\left({{0.2755731922_3985890653 E - 6}\ {{lambda}^{10}}\ {{e}^{-{{1.0}\  lambda}}}}\right)}}}
(2)
Type: Expression(Float)
fricas
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 BOOT::|*2;optimFun;1;initial| is undefined.

As an aside: can we calculate an elasticity from logistic regression?

fricas
Y:=exp(alpha+betak*Pk)/(1+exp(alpha+betak*Pk))

\label{eq3}{{e}^{{Pk \  betak}+ alpha}}\over{{{e}^{{Pk \  betak}+ alpha}}+ 1}(3)
Type: Expression(Integer)
fricas
D(Y,Pk)/Y

\label{eq4}betak \over{{{e}^{{Pk \  betak}+ alpha}}+ 1}(4)
Type: Expression(Integer)
fricas
D(Y,Pk)/Y/(1-Y)

\label{eq5}betak(5)
Type: Expression(Integer)




  Subject:   Be Bold !!
  ( 14 subscribers )  
Please rate this page: