Re: [Axiom-math] (no subject)

[Martin Rubey]

C. Frangos writes:
 > 
 > 17 June 2005
 > 
 > Some time ago, I managed to get axiom running under suse linux 7.2 (graphics
 > not working). Recently I downloaded the latest version and got it running
 > under suse linux 9.2 (with working graphics) on another machine.

great!

 > I am trying to convert some of my Matlab m-files to axiom. 

even better!

 > Any assistance with the following will be much appreciated:

let's try!

 > (1) Is there a short way of declaring all variables in an axiom function to
 >     be local ?


If you compile your files, this is automatically the case. However, it seems
that you want to use .input files, not .spad files?

Thinking about it, it should be automatically the case in .input files,
too. Could you send an example?

 >     In Matlab this is automatically the case and global variables have
 >     to be specifically declared. 

Mostly the same in Axiom. See section 6.16 of the Axiom Book.
 
 >     In axiom it seems that its the opposite, creating the possibility
 >     of errors. Is this a design error in an otherwise very
 >     advanced language ?
 
 > (2) In Matlab I have a main function (in an m-file) calling another
 >     function, which in turn calls another function etc. Its not clear
 >     from the axiom book how to implement this. Does one have to use
 >     )read fun1.input, )read fun2.input, ...etc, and then fun1(); ?

No, functions can call each other within on input file.

 >     Any papers, reports or tutorials on the details of axiom programming
 >     involving multiple functions defined in separate files, calling
 >     each other, etc as described above would be helpful.

See below. I'll include my latest input file that deals with garch processes.
 
 > (3) Is there a setting in axiom to avoid having to use underscore _,
 >     to insert blank lines between program lines in files defining
 >     functions ? This seems awkward and error prone.

You don't have to do this usually. Only if you want to continue a line.


 > (4) I am using the Matlab symbolic toolbox (Maple engine). Its not
 >     clear from the axiom book how to typecast a unit matrix (integer
 >     or float) to type expression or variable in order to add it to
 >     another matrix containing expressions in the variable/symbol x,
 >     for example: A:= I + matrixins;

Don't do this. You want to use the + from MATRIX INT, for example.

m1 := [[1,2],[3,4]]
m2 := [[x,y],[u,v]]

m1+m2 should add the two without complaining

 > (5) System commands like )read, etc don't seem to work within
 >     functions ?

No.

 > (6) Is there an axiom equivalent to Matlab's pause command ?

Don't know.
 
 > (7) I tried using a list to return multiple variables from a
 > function. However, this turns out to be of type any. These variables can be
 > displayed but not used in operations any further. This means that you have
 > to remember the types of each variable in the list.  Then, manually in the
 > calling function extract and explicitly typecast each variable in the
 > returned list. Is there an easier way of doing this ??

They shouldn't be converted to type ANY. Often it helps to declare your
functions. See below.

 > (8) Does the current axiom distribution with working graphics have to be
 > recompiled on my older suse linux 7.2 machine (first one took about 12
 > hours) or can I somehow copy the directory of the latest axiom via ethernet
 > from my newer suse linux 9.2 machine ?

Yes, that should work. (If the paths are the same...)

Here comes a .input from my research.

epsilon contains a garch process after you said

)re garch.input

from within axiom


Hope this helps a little bit. I'll try to detail (2) coming monday.

Martin

-------------------------------------------------------------------------------


-- n:=100; max:=1.2*Lik2(n,w,a,b,eps); 
draw((w:DoubleFloat,a:DoubleFloat):DoubleFloat +-> 
min(Lik2(n,w,a,b,eps),max),0..3,0..4)
-- draw((w:DoubleFloat):DoubleFloat +-> Lik2(n,w,a,b,eps),0..0.5)
-- draw([[i,eps.i]::Point Float for i in 1..1000])


-- eta:=normal(0,1)$RFDIST

eta:=uniform(-sqrt(3),sqrt(3))$RFDIST

-- eta aus -1, 1, 3, 5 -> eps > w

-- eta():INT == 2*(uniform(0..3)$RIDIST)()-1 

w := 0.3356; a := 2; b := 0;

-- without this line it's not a Garch process!
)set fun cache 2 epsilon 
 
epsilon(0) == 0;
epsilon(t:INT):FLOAT == sigma(t)*eta();

)set fun cache 2 sigma

sigma(0) == 0;
sigma(t:INT):FLOAT == 
  sqrt(w+a*epsilon(t-1)^2 + b*sigma(t-1)^2)

precision(20);
eps := [epsilon(t) for t in 1..1000];

)set fun cache 2 sigma2

sigma2(0,w,a,b,eps) == 0;
sigma2(t:INT,w:FLOAT,a:FLOAT,b:FLOAT,eps:LIST FLOAT):FLOAT == 
  w + a*(eps.t)^2 + b*sigma2(t-1,w,a,b,eps)

sigma2p(0,w,x,y,eps) == 0;
sigma2p(t:INT,w:FLOAT,x:FLOAT,y:FLOAT,eps:LIST FLOAT):FLOAT == 
  w*(1 + x*(eps.t)^2 + y*sigma2p(t-1,w,x,y,eps))

Likp(n:INT,x:FLOAT,eps:LIST FLOAT):FLOAT ==
  log(1/n*reduce(+,[(eps.(t+1))^2/(1+x*(eps.t)^2) for t in 1..n])) _ 
   -1/n*reduce(+,[log(1/(1+x*(eps.t)^2)) for t in 1..n])+1

lik(t:INT,w:FLOAT,a:FLOAT,b:FLOAT,eps:LIST FLOAT):Union(FLOAT, "fail") ==
  s := sigma2(t,w,a,b,eps)
  if s <= 0
  then "fail"
  else log(s)+(eps.(t+1))^2/s

lik2(t:INT,w:FLOAT,a:FLOAT,b:FLOAT,eps:LIST FLOAT):FLOAT == 
  (l := lik(t,w,a,b,eps)) case "fail" => 1000000.0
  l

lik3(t,w,a,eps) ==
  s := w+a*eps.t^2
  if s <= 0
  then "fail"
  else log(s)+(eps.(t+1))^2/s

Lik(n:INT,w:FLOAT,a:FLOAT,b:FLOAT,eps:LIST FLOAT):Union(FLOAT, "fail") == 
  ll := 0.0
  for t in 1..n repeat
    (l := lik(t,w,a,b,eps)) case "fail" => return "fail"
    ll := ll + l::Float
  ll/n

Lik2(n:INT,w:FLOAT,a:FLOAT,b:FLOAT,eps:LIST FLOAT):FLOAT == 
  (ll := Lik(n,w,a,b,eps)) case "fail" => 1000000.0
  ll

Lik3(n,w,a,eps) == 
  1/n*reduce(+,[lik3(t,w,a,eps) for t in 1..n])

minInd(l: List List Float): List List INT == 
  m := l.(1,1)
  ind := [[1,1]]
  for i in 1..#l repeat
    for j in 1..#l.i repeat
      if l.i.j = m then 
        ind := cons([i,j], ind)
      if l.i.j < m then 
        m := l.i.j
        ind := [[i,j]]
  ind

minInd1(l: List Float): List INT == 
  m := l.1
  ind := [1]
  for i in 1..#l repeat
    if l.i = m then 
      ind := cons(i, ind)
    if l.i < m then 
      m := l.i
      ind := [i]
  ind

Likp2(n,x,eps) ==
  log(1/n*reduce(+,[(eps.(t+1))^2/(1+x*(eps.t)^2) for t in 1..n])) _ 
   -1/n*reduce(+,[log(1/(1+x*(eps.t)^2)) for t in 1..n])+1

Likp2a(n,x,eps) == log(1/n*reduce(+,[(eps.(t+1))^2/(1+x*(eps.t)^2) for t in 
1..n]))

Likp2b(n,x,eps) == -1/n*reduce(+,[log(1/(1+x*(eps.t)^2)) for t in 1..n])+1

LikpLim(n,eps) == LikpLima(n,eps) + LikpLimb(n,eps) + 1

LikpLima(n,eps) == 
  log(reduce(+,[eps.(i+1)^2/eps.i^2 for i in 1..n])/n)

LikpLimb(n,eps) == 
  log(reduce(*,[eps.i^2 for i in 1..n]))/n

e:=[concat("e",i::String)::Symbol::EXPR INT for i in 1..20];

f(2)==e.1^4*e.3^2+e.2^6
f(n)==(n-1)*(reduce(*,[e.i for i in 1..n-1])^4*e.(n+1)^2+e.n^4*f(n-1)/(n-2))

g(2)==e.2^2+e.3^2
g(n)==g(n-1)*e.(n+1)^2+reduce(*,[e.i^2 for i in 1..n-2])*e.n^2*e.(n+1)^2

-- n := 5; factor(reduce(+, select(e +-> leadingCoefficient (e) = n-1, 
monomials leadingCoefficient(numer(D(Likp2(n,x,[e1,e2,e3,e4,e5,e6]),x)::FRAC 
POLY INT::FRAC UP(x, POLY INT)))))::DMP([e1,e2,e3,e4,e5,e6], INT))-f(n)

-- n := 5; factor(reduce(+, select(e +-> leadingCoefficient (e) = -1, monomials 
leadingCoefficient(numer(D(Likp2(n,x,[e1,e2,e3,e4,e5,e6]),x)::FRAC POLY 
INT::FRAC UP(x, POLY INT)))))::DMP([e1,e2,e3,e4,e5,e6], INT))



-- n:=normal(0,1)$RFDIST
-- 
-- epsilon(w,a,b,0) == 0
-- epsilon(w:FLOAT,a:FLOAT,b:FLOAT,t:INT):FLOAT == sigma(w,a,b,t)*rnormal()
-- 
-- sigma(w,a,b,0) == 0
-- sigma(w:FLOAT,a:FLOAT,b:FLOAT,t:INT):FLOAT == 
sqrt(w+a*epsilon(w,a,b,t-1)^2+b*sigma(w,a,b,t-1)^2) 
-- 
-- sigma2(0,w,a,b,eps) == 0
-- sigma2(t:INT,w:FLOAT,a:FLOAT,b:FLOAT,eps:LIST FLOAT):FLOAT == 
w+a*(eps.t)^2+b*sigma2(t-1,w,a,b,eps)^2
-- 
-- lik(t:INT,w:FLOAT,a:FLOAT,b:FLOAT,eps:LIST FLOAT):FLOAT == 
log(sigma2(t,w,a,b,eps))+(eps.(t+1))^2/sigma2(t,w,a,b,eps)
-- Lik(n:INT,w:FLOAT,a:FLOAT,b:FLOAT,eps:LIST FLOAT):FLOAT == 
reduce(+,[lik(t,w,a,b,eps) for t in 1..n])/n
-- 
-- precision(10)
-- eps := [epsilon(1,1/2,1/2,t) for t in 1..10]