From: Pascal J. Bourguignon
Subject: Re: Lisp and Symbolic Integration
Date: 
Message-ID: <7cwslvzsk2.fsf@pbourguignon.anevia.com>
Francogrex <······@grex.org> writes:

> Dear Lispers,
> First introduction: I am new to Lisp programming (Common Lisp), I
> discovered it about a week while looking for a language that
> especially supports well symbolic mathematical analysis. 

Welcome! 

> I have used C+
> + to write programs before but never really warmed to that language. 

Who can?


> I
> also know how to program well with S+ and R two statistical computing
> programs (but "specific" languages), said to be based loosely on Lisp.
> I immediately loved Common Lisp, despite the "strange" syntax at
> first. I tried some implementations but decided that the one I'm most
> confortable with is GCL, it's nice and easy to use (I also tried CLISP
> but was put off by floating point error messages reminiscent of the C
> languages, so I dropped it). 


It may be not error messages, but warnings, about underflow?  There's
an option to silent. Have a look 
at: CUSTOM:*WARN-ON-FLOATING-POINT-CONTAGION*
on: http://clisp.cons.org/impnotes/num-concepts.html

You may want to try again clisp, knowing that you can set the length
of the mantissa of long-floats (and then don't have any limit on the
exponent too):

C/USER[44]> (SETF (EXT:LONG-FLOAT-DIGITS) 1000)
1000
C/USER[45]> (length (prin1-to-string pi))
310
C/USER[46]> (SETF (EXT:LONG-FLOAT-DIGITS) 4000)
4000
C/USER[47]> (length (prin1-to-string pi))
1208
C/USER[48]> 



> I also didn't like emacs much so I
> settled for my own editors which has extra features to match
> parenthesis and highlight keywords etc (called zion). I have already
> converted several of my C++ and R programs into lisp and they look
> really neat. 

When I converted my programs from C++, the lisp programs were between
1/10 and 1/30 of size of the equivalent C++ programs.  Did you observe
the same for yours?


> I will ask one question but I'm not sure it's something
> available: Does any of you guys know where to find a good webstore of
> Lisp programs? I'm particularly looking for the "Hermite Reduction"
> algorithm for symbolic integration, and wonder if someone has already
> written a lisp code for it. Thanks.

Have a look at:
http://www.cliki.net/
http://common-lisp.net/
http://www.cl-user.net/asp/root-dir

For mathematical functions, I'd guess the sources of maxima
would be a good repository.

http://www.cliki.net/Maxima
http://www.cliki.net/cl-mathstats
http://www.cl-user.net/asp/search?search=math


-- 
__Pascal Bourguignon__

From: Robert Dodier
Subject: Re: Lisp and Symbolic Integration
Date: 
Message-ID: <1d49664c-9d13-4332-af84-71a6d6a81e4d@34g2000hsh.googlegroups.com>
Francogrex wrote:

> Hi, thanks. Actually it's a floating-point-overflow errors. Funny that
> I don't get those while using GCL and my results in GCL are the same
> as those in R/S+ (meaning GCL does the good job).

A close reading of the Common Lisp Hyperspec seems to indicate [1]
that indefinite or infinite floats do not exist in CL.
GCL is just being sloppy: it is not checking the result of a
floating point operation to ensure that it is definite and finite.
Some other Lisp implementations, such as Clisp, interpret the spec
narrowly and banish indefinite or infinite floats; still others allow
the
programmer to permit or banish such floats, according to some flags.

> I really get the feeling that GCL is more performant and easier to handle.

GCL has various idiosyncrasies, and development seems sporadic.
I'll recommend SBCL instead. It is a more careful implementation of
the CL spec, and it is almost as fast as GCL. You can set some flags
to make SBCL allow floating point not-a-number and infinity.

FWIW

Robert Dodier

[1] Seems, that is, to those who are sufficiently enlightened;
one achieves enlightenment by reading the spec many times and
meditating deeply on its mysteries.
From: Robert Dodier
Subject: Re: Lisp and Symbolic Integration
Date: 
Message-ID: <d5d32ad2-e889-4455-b439-56ca0e9ea467@p25g2000hsf.googlegroups.com>
On May 15, 9:32 am, Francogrex <······@grex.org> wrote:

> Maybe, but I don't think GCL is being sloppy, because R and S-plus
> work the same on this one as GCL and output exactly the same results.
> If GCL is sloppy then so is R (which seems strange).

That's beside the point; R is not an implementation of
Common Lisp, so it is not required to obey the same
rules about floating point numbers.

Robert Dodier
From: Lars Rune Nøstdal
Subject: Re: Lisp and Symbolic Integration
Date: 
Message-ID: <482c969c$0$2323$c83e3ef6@nn1-read.tele2.net>
Francogrex wrote:
> On May 15, 7:53 pm, Robert Dodier <·············@gmail.com> wrote:
>> On May 15, 9:32 am, Francogrex <······@grex.org> wrote:
>>> Maybe, but I don't think GCL is being sloppy, because R and S-plus
>>> work the same on this one as GCL and output exactly the same results.
>>> If GCL is sloppy then so is R (which seems strange).
>> That's beside the point; R is not an implementation of
>> Common Lisp, so it is not required to obey the same
>> rules about floating point numbers.
> 
> Listen guys, you are the experts since you've been using lisp longer
> than I have. But compliant or not (though my version 2.6.7 is ansi
> compliant), GCL does what I need; when I type (exp 700) in GCL or
> Exp[700]//N in Mathematica I get what I need and I need large numbers
> and divisions by large numbers. Clisp isn't much help there for me.
> Just so that I know, is GCL looked down upon here?

does it matter? (ok, yeah .. i guess sometimes it does if you want
help in the future)

fwiw. i'd stick with sbcl or clisp or something one could be sure
worked in the future

i mean; if gcl "suddenly" becomes more compliant wrt. the Common
Lisp spec it would mean that the software you've written will
break .. wouldn't it?

i think i'd look for, perhaps, a lisp library that does what
you need regarding numbers

-- 
Lars Rune N�stdal
http://nostdal.org/
From: John Thingstad
Subject: Re: Lisp and Symbolic Integration
Date: 
Message-ID: <op.ua7ozficut4oq5@pandora.alfanett.no>
P� Thu, 15 May 2008 22:01:29 +0200, skrev Lars Rune N�stdal  
<···········@gmail.com>:

>
> i mean; if gcl "suddenly" becomes more compliant wrt. the Common
> Lisp spec it would mean that the software you've written will
> break .. wouldn't it?

Naw, it means it is compliant with CLTL2 by Steele.
It is unlikely to break because it becomes ANSI compliant.

That said CLisp uses the GNU multi-precision library too (like  
Mathematica) snd has arbitrary precision float.

--------------
John Thingstad
From: Barry Fishman
Subject: Re: Lisp and Symbolic Integration
Date: 
Message-ID: <m3wslvi578.fsf@barry_fishman.acm.org>
Francogrex <······@grex.org> writes:

> On May 15, 7:53 pm, Robert Dodier <·············@gmail.com> wrote:
>> On May 15, 9:32 am, Francogrex <······@grex.org> wrote:
>> > Maybe, but I don't think GCL is being sloppy, because R and S-plus
>> > work the same on this one as GCL and output exactly the same results.
>> > If GCL is sloppy then so is R (which seems strange).
>> That's beside the point; R is not an implementation of
>> Common Lisp, so it is not required to obey the same
>> rules about floating point numbers.
>
> Listen guys, you are the experts since you've been using lisp longer
> than I have. But compliant or not (though my version 2.6.7 is ansi
> compliant), GCL does what I need; when I type (exp 700) in GCL or
> Exp[700]//N in Mathematica I get what I need and I need large numbers
> and divisions by large numbers. Clisp isn't much help there for me.
> Just so that I know, is GCL looked down upon here?

As far as large numbers are concerned:

$ gcl
GCL (GNU Common Lisp)  2.7.0 ANSI    Jan 11 2008 19:03:17
Source License: LGPL(gcl,gmp,pargcl), GPL(unexec,bfd,xgcl)
Binary License:  GPL due to GPL'ed components: (XGCL READLINE BFD UNEXEC)
Modifications of this banner must retain notice of a compatible license
Dedicated to the memory of W. Schelter

Use (help) to get some basic information on how to use GCL.

Temporary directory for compiler files set to /tmp/

>most-positive-short-float

3.40282347S38

>most-positive-single-float

1.7976931348623157E308

>most-positive-double-float

1.7976931348623157E308

>most-positive-long-float

1.7976931348623157E308

>(type-of (exp 400))

LONG-FLOAT

>>(exp 800l0)

Error:
Fast links are on: do (si::use-fast-links nil) for debugging
Signalled by EVAL.
Condition in EVAL [or a callee]: INTERNAL-SIMPLE-FLOATING-POINT-OVERFLOW: Arithmetic error when performing EXP on 800.0:

Broken at EVAL.  Type :H for Help.
 1 (Continue) Return to top level.
COMMON-LISP-USER>>

So GCL is using double precision for its default single-precision
floating point type, while most other lisps using 32bit single precision.

However:

$ clisp
  i i i i i i i       ooooo    o        ooooooo   ooooo   ooooo
  I I I I I I I      8     8   8           8     8     o  8    8
  I  \ `+' /  I      8         8           8     8        8    8
   \  `-+-'  /       8         8           8      ooooo   8oooo
    `-__|__-'        8         8           8           8  8
        |            8     o   8           8     o     8  8
  ------+------       ooooo    8oooooo  ooo8ooo   ooooo   8

Welcome to GNU CLISP 2.45 (2008-04-04) <http://clisp.cons.org/>

Copyright (c) Bruno Haible, Michael Stoll 1992, 1993
Copyright (c) Bruno Haible, Marcus Daniels 1994-1997
Copyright (c) Bruno Haible, Pierpaolo Bernardi, Sam Steingold 1998
Copyright (c) Bruno Haible, Sam Steingold 1999-2000
Copyright (c) Sam Steingold, Bruno Haible 2001-2008

Type :h and hit Enter for context help.

;; Loading file /home/barry/.clisprc.lisp ...
;;  Loading file /home/util/lisp/site/asdf.fas ...
;;  Loaded file /home/util/lisp/site/asdf.fas
;;  Loading file /home/util/lisp/site/defsystem.fas ...
;;  Loaded file /home/util/lisp/site/defsystem.fas
;; Loaded file /home/barry/.clisprc.lisp
[1]> most-positive-short-float
3.4028s38
[2]> most-positive-single-float
3.4028235E38
[3]> most-positive-double-float
1.7976931348623157d308
[4]> most-positive-long-float
8.8080652584198167656L646456992
[5]> (exp 400d0)
5.221469689764144d173
[6]> (exp 800l0)
2.7263745721125665673L347

So CLISP can actually handle larger numbers if they are expressed in
long-float format.

If you want to keep you Lisp options open, you just need to avoid
default conversions when needing large numbers.

-- 
Barry Fishman
From: Pascal J. Bourguignon
Subject: Re: Lisp and Symbolic Integration
Date: 
Message-ID: <7cprrmxwfq.fsf@pbourguignon.anevia.com>
Francogrex <······@grex.org> writes:

> On May 16, 12:16 am, Barry Fishman <·············@acm.org> wrote:
>> So CLISP can actually handle larger numbers if they are expressed in
>> long-float format.
>> If you want to keep you Lisp options open, you just need to avoid
>> default conversions when needing large numbers.
>
> Thanks, now I see why this happens. Can we globally set all double to
> long floats in clisp

You can set *read-default-float-format* so that when you read a
floating point number that has no exponent specifying the type of
floating point, it is read as the type you want.

(setf *read-default-float-format* 'long-float)

Note that by default, it's simple-float, not even double-float.

(It may happen that simple-float = double-float = 64-bit IEEE754
in some implementations, thought).

> so that it behaves like GCL? 

In clisp, with the default  *read-default-float-format*, 

C/USER[87]> (type-of 1.2)
SINGLE-FLOAT
C/USER[88]> (type-of 1.2e2)
SINGLE-FLOAT
C/USER[89]> (type-of 1.2d2)
DOUBLE-FLOAT
C/USER[90]> (type-of 1.2l2)
LONG-FLOAT
C/USER[91]> 

so perhaps what you're complaining is that you're using 32-bit floats
in clisp, vs. 64-bit floats in gcl, by default?




> A little bit more practical, is there
> anybody who can help me fix
> those double/long float issues in my program below so that it would be
> possible to run it on clisp?
> I have commented the program and also did some minor formatting.
>
> ;; The while macro
> (defmacro while (condition &rest body)
>         (let ((var (gensym)))
>             `(do ((,var nil (progn ,@body)))
>                  ((null ,condition) ,var))))

(defmacro while (condition &body body)
   `(do () ((not ,condition)) ,@body))


> ;;The factorial function
> (defun factorial (x) (if (eql x 0) 1 (* x (factorial (- x 1)))))

Use = to compare numbers.  (eql 0 0.0) --> NIL
Better, in the case of factorial use <=

(defun factorial (x)
  (if (<= x 0.0) 1.0 (* x (factorial (- x 1.0)))))

> ;;The poisson PDF
> (defun dpois(N L)
> 		(mapcar #'(lambda (x) (/( *(expt L x)  (exp (- L))) (factorial x)))
> N))
> ;;The sum function
> (defun sum(zazlist) (loop for x in zazlist sum x))
> ;;The mean function

> (defun mean(zazlist)
> 					(setq sum (loop for x in zazlist sum x))
> 					(setq m (/ sum (length zazlist))))

Use local variables, not undefined things.  

(defun mean (zazlist)
   (let ((sum  (loop for x in zazlist sum x))
         (mean (/ sum (length zazlist))))
      mean))


Use your own abstractions!

(defun mean (zazlist)
   (/ (sum zazlist) (length zazlist)))

       
> ;;Set eta (the tuning parameter)
> (setq eta (list 1 1 1))
> ;;Set the counter.
> (setq count 0)

In a program, you should use DEFPARAMETER (or DEFVAR) to define
variables.  There's no telling what SETQ will do, it can change from
one implementation to the other. Also, DEFPARAMETER/DEFVAR allow you
to document your variables:

(defparameter *eta* (list 1 1 1) "The tuning parameter.")

So (describe '*eta*) now prints:

------------------------------------------------------------------------
*ETA* is the symbol *ETA*, lies in #<PACKAGE COMMON-LISP-USER>, is
accessible in 1 package COMMON-LISP-USER, a variable declared
SPECIAL, value: (1 1 1), has 1 property SYSTEM::DOC.

Documentation as a VARIABLE:
The tuning parameter.
For more information, evaluate (SYMBOL-PLIST '*ETA*).

 #<PACKAGE COMMON-LISP-USER> is the package named COMMON-LISP-USER. It has 2 nicknames CL-USER, USER.
 It imports the external symbols of 2 packages COM.INFORMATIMAGO.PJB,
 COMMON-LISP and exports 7 symbols, but no package uses these exports.

 (1 1 1) is a list of length 3.

Documentation:
VARIABLE:
"The tuning parameter."
------------------------------------------------------------------------


> ;;The Expectation-Maximization (EM) algorithm function to find the two
> means of the poisson mixtures (p2 & p3)
> ;;and the mixing probability p1.

By the way, you should also better put the documentation of your
function as documentation string rather than ;-comments.  These
;-comments are just discarded at read-time, they're totally useless.

(defun sum (zazlist)
    "
Return:     Σ x
         x ∈ zazlist
"
    (loop :for x :in zazlist :sum x))


So now: (progn (princ (documentation 'sum 'function)) (values))
will print:
------------------------------------------------------------------------
Return:     Σ x
         x ∈ zazlist
------------------------------------------------------------------------


> (defun emalgo(N  p1 p2 p3)
> 			(while (or (> (nth 0 eta) 0.01 ) (> (nth 1 eta) 0.01 ) (> (nth 2
> eta) 0.01 ))
> 					(incf count)
> 					(setq op1 p1)
> 					(setq op2 p2)
> 					(setq op3 p3)
> 					;;The Expectation step
> 					(setq EZ (mapcar #'/ (mapcar #'(lambda(x) (* p1 x)) (dpois N
> p2))
> 									(mapcar #'+ ( mapcar #'(lambda(x) (* p1 x)) (dpois N p2))
> (mapcar
> 									#'(lambda(x) (* (- 1 p1) x)) (dpois N p3)))))
> 					;;The Maximization steps
> 					(setq p1 (mean EZ))
> 					(setq p2 (/(sum (mapcar #'* EZ N )) (sum EZ)))
> 					(setq p3 (/(sum (mapcar #'* (mapcar '(lambda(x) (- 1 x)) EZ) N ))
> (sum
> 					(mapcar '(lambda(x) (- 1 x)) EZ))))
> 					(setq eta  (mapcar #'abs (mapcar #'- (list op1 op2 op3) (list p1
> p2 p3))))
> ;;Print out eta and the intermediate results of the p1,p2,p3
> estimates.
> (print (list "eta:" eta))
> (print (list "count:" count "p" (list p1 p2 p3)))))

No problem in using PRINT, while debugging, notably in combination
with backquote and comma:

(print `(eta= ,eta))

(let ((p (list 1.2 2.3 3.4))
      (count 42))
  (print `(count= ,count p= ,@p)))
------------------------------------------------------------------------
(COUNT= 42 P= 1.2 2.3 3.4) 
------------------------------------------------------------------------

but for nicer output, you could use FORMAT:

(let ((p (list 1.2 2.3 3.4))
      (count 42))
  (format t "~&Count: ~8D~%P:     ~{~8F~^,~}~%" count p))
------------------------------------------------------------------------
Count:       42
P:          1.2,     2.3,     3.4
------------------------------------------------------------------------


> ;;Run an example:
> (setq N (list 23 24 78 65 34 98 12 23 54 34 ))
> (emalgo N 0.2 20 10)

Otherwise, as long as the arguments are integer, the results of most
arithmetic operations will be integer, or rational.  Not floating
point numbers.

If you are sure you want floating point results, it might be
interesting to start to use floating points earlier in computation,
before the numbers become bignums, or rational (which may become ratio
of bignums rather quickly).



(setf *read-default-float-format* 'double-float)

(defun factorial (x)
  "Factorial in floating point numbers."
  (declare (type #.*read-default-float-format* x))
  (if (<= x 0.0) 1.0 (* x (factorial (- x 1.0)))))

(defun dpois(N L)
  "The poisson PDF"
  (mapcar (lambda (x) (/ (* (expt L x) (exp (- L))) (factorial x)))
          N))

(defun sum (zazlist)
    "
Return:     Σ x
         x ∈ zazlist
"
    (loop :for x :in zazlist :sum x))

(defun mean (zazlist)
   (/ (sum zazlist) (length zazlist)))


;; In general, it's better to avoid global variables.
;; Often you can gather the general parameters of a problem into an
;; object or a structure:

(defstruct (problem-parameters (:conc-name problem-))
  (eta   (list 1.0 1.0 1.0))
  (count 0))



;; For most numerical algorithms, I'd rather use vectors than lists.
;; There is no list in lisp. Yes, it's paradoxical.  You could  still
;; invent a list-of type such as:

(deftype list-of (element-type)
  (let ((predicate-name (gensym)))
    (compile predicate-name
             `(lambda (object)
                (and (listp object)                      ; a list
                     (list-length object) ; a proper list
                     (every (lambda (element) (typep element ',element-type))
                            object))))
    `(satisfies ,predicate-name)))


;; C/USER[114]> (typep '(1 2 3) '(list-of integer))
;; T
;; C/USER[115]> (typep '(1.0 2 3) '(list-of integer))
;; NIL
;; C/USER[116]> (typep '("1.0" "2" "3") '(list-of string))
;; T
;; C/USER[117]> (typep '("1.0" 2 "3") '(list-of string))
;; NIL

;; But most compilers wont handle it with the same efficacity nor
;; expediency than a vector type such as (vector double-float).


(defun emalgo (problem  N  p1 p2 p3)
  (declare (type problem-parameters problem)
           (type #.*read-default-float-format* p1 p2 p3)
           (type (list-of #.*read-default-float-format*) N))
  (flet ((1-x (x) (- 1.0 x)))
    (let ((eta (problem-eta problem)))
      (loop :while (some (lambda (x)  (< 0.01 x)) eta) :do
         (incf (problem-count problem))
         (let ((op1 p1)
               (op2 p2)
               (op3 p3)
               ;;The Expectation step:
               (EZ (mapcar (function /)
                           (mapcar (lambda(x) (* p1 x))
                                   (dpois N p2))
                           (mapcar (function +)
                                   (mapcar (lambda(x) (* p1 x))
                                           (dpois N p2))
                                   (mapcar (lambda(x) (* (- 1.0 p1) x))
                                           (dpois N p3))))))
        
           ;;The Maximization steps
           (setf p1 (mean EZ))
           (setf p2 (/ (sum (mapcar (function *) EZ N )) (sum EZ)))
           (setf p3 (/ (sum (mapcar (function *) (mapcar (function 1-x) EZ) N))
                       (sum (mapcar (function 1-x) EZ))))
           (setf eta (mapcar (function abs)
                             (mapcar (function -)
                                     (list op1 op2 op3)
                                     (list  p1  p2  p3))))))
      (setf (problem-eta problem) eta))))



C/USER[137]> (let ((problem  (make-problem-parameters)))
               (prog1 (emalgo problem
                              (list 23 24 78 65 34 98 12 23 54 34 )
                              0.2 20.0 10.0)
                 (print problem)))
Prints:
------------------------------------------------------------------------
#S(PROBLEM-PARAMETERS :ETA (2.41653767321659E-5 0.0019710634938405747 6.493578587161153E-4) :COUNT 5) 
------------------------------------------------------------------------

--> (2.41653767321659E-5 0.0019710634938405747 6.493578587161153E-4)


If double-float are not good enough, you will have to change the
setting of the *read-default-float-format* variable, and
reload/recompile the source file.

When using #+clisp long-float, remember to set their length you want.

-- 
__Pascal Bourguignon__
From: Barry Fishman
Subject: Re: Lisp and Symbolic Integration
Date: 
Message-ID: <m3hccyfbnk.fsf@barry_fishman.acm.org>
···@informatimago.com (Pascal J. Bourguignon) writes:
> In a program, you should use DEFPARAMETER (or DEFVAR) to define
> variables.  There's no telling what SETQ will do, it can change from
> one implementation to the other. Also, DEFPARAMETER/DEFVAR allow you
> to document your variables:
>
> (defparameter *eta* (list 1 1 1) "The tuning parameter.")

Eeek. I didn't have that of my response.  Global values always should
begin and end with "*".  Besides making their special behavior more
evident, they stand out better in your code.  This also avoids the
problem with COUNT as a variable, it should really be *COUNT* anyway.
The values still need to be given in floating point format.

(defparameter *eta* (list 1.0 1.0 1.0) "The tuning parameter.")

-- 
Barry Fishman
From: Barry Fishman
Subject: Re: Lisp and Symbolic Integration
Date: 
Message-ID: <m3mymqs2bf.fsf@barry_fishman.acm.org>
--=-=-=
Content-Type: text/plain; charset=utf-8
Content-Transfer-Encoding: 8bit

Francogrex <······@grex.org> writes:

> On May 16, 12:16 am, Barry Fishman <·············@acm.org> wrote:
>> So CLISP can actually handle larger numbers if they are expressed in
>> long-float format.  If you want to keep you Lisp options open, you
>> just need to avoid default conversions when needing large numbers.
>
> Thanks, now I see why this happens. Can we globally set all double to
> long floats in clisp so that it behaves like GCL?

You can do:

(setf *read-default-float-format* 'long-float)

This will make the float type for literal floating point numbers (like
1.0 or 1.0e3) in your program LONG-FLOAT.  It will not make "(exp 400)"
work because the 400 will be considered an integer argument to exp and
it will still end up coerced to single precision.  Therefore you must be
careful to make sure numbers that are not always to be considered
integers are entered in floating point format.  In general expresions
involving long-float numbers will result in long-float values, although
there might be occasions you need to force the conversion with something
like:

(let ((x 400))
   (exp (coerce x 'long-float)))

>                                                      A little bit more
> practical, is there anybody who can help me fix those double/long
> float issues in my program below so that it would be possible to run
> it on clisp?  I have commented the program and also did some minor
> formatting.
>
...

> ;;The factorial function
> (defun factorial (x) (if (eql x 0) 1 (* x (factorial (- x 1)))))

Although this will work with integer arguments, EQL is not the best
choice for numeric comparisons.  It would be better to use "(= x 0)"
or even in this case "(zerop x)".

Also multiple recursive calls can cause problems, so LOOP style
programming is sometimes needed.  In your example code this isn't an
issue.

>...
> ;;Set eta (the tuning parameter)
> (setq eta (list 1 1 1))

These are treated as floating point numbers, so they should be
entered that way.

(setq eta (list 1.0 1.0 1.0))

Of if you do not set *read-default-float-format*:

(setq eta (list 1.0l0 1.0l0 1.0l0))

> (mapcar #'+ ( mapcar #'(lambda(x) (* p1 x)) (dpois N p2))
>   (mapcar #'(lambda(x) (* (- 1 p1) x)) (dpois N p3)))))

You can leave the #' prefix off of (lambda (...) ...)

> ;;Run an example:
> (setq N (list 23 24 78 65 34 98 12 23 54 34 ))
> (emalgo N 0.2 20 10)

N is treated as a list of integers so it need not be changed, but,
I think it is necessary to change the arguments to emalgo:

(emalgo N 0.2 20.0 10.0)

again assuming that *read-default-float-format* is set.

I have attached my own quick reformatted and changed code. (I made
it into a unix style clisp script but you can just remove the initial #!
line to LOAD it normally.)  It seems to give the same results as your
program but with a few more digits displayed.


--=-=-=
Content-Type: x-lisp
Content-Disposition: attachment; filename=tryme.lisp
Content-Description: emalgo script

#! /usr/bin/env clisp
(setf *read-default-float-format* 'long-float)

;; The while macro
(defmacro while (condition &rest body)
  (let ((var (gensym)))
    `(do ((,var nil (progn ,@body)))
	 ((null ,condition) ,var))))

;;The factorial function
(defun factorial (x) (if (zerop x) 1 (* x (factorial (- x 1)))))

;;The poisson PDF
(defun dpois(N L)
  (mapcar (lambda (x) (/ (* (expt L x) (exp (- L))) (factorial x)))
	  N))

;;The sum function
(defun sum (zazlist) (loop for x in zazlist sum x))

;;The mean function
(defun mean(zazlist)
  (setq sum (loop for x in zazlist sum x))
  (setq m (/ sum (length zazlist))))

;;Set eta (the tuning parameter)
(setq eta (list 1.0 1.0 1.0))

;;Set the counter.
(setq count 0)

;; The Expectation-Maximization (EM) algorithm function to find the
;; two means of the poisson mixtures (p2 & p3) and the mixing
;; probability p1.
(defun emalgo(N  p1 p2 p3)
  (while (or (> (nth 0 eta) 0.01)
	     (> (nth 1 eta) 0.01)
	     (> (nth 2 eta) 0.01))
    (incf count)
    (setq op1 p1)
    (setq op2 p2)
    (setq op3 p3)
    ;;The Expectation step
    (setq EZ (mapcar #'/ (mapcar (lambda (x) (* p1 x))
				 (dpois N p2))
		     (mapcar #'+ (mapcar (lambda (x) (* p1 x))
					 (dpois N p2))
			     (mapcar (lambda (x) (* (- 1 p1) x))
				     (dpois N p3)))))

    ;;The Maximization steps
    (setq p1 (mean EZ))
    (setq p2 (/ (sum (mapcar #'* EZ N )) (sum EZ)))
    (setq p3 (/ (sum (mapcar #'* (mapcar (lambda (x) (- 1.0 x)) EZ) N))
	        (sum (mapcar (lambda (x) (- 1.0 x))
			     EZ))))
    (setq eta (mapcar #'abs (mapcar #'- (list op1 op2 op3)
				    (list p1 p2 p3))))

    ;;Print out eta and the intermediate results of the p1,p2,p3 estimates.
    (print (list "eta:" eta))
    (print (list "count:" count "p" (list p1 p2 p3)))))

;;Run an example:
(setq N (list 23 24 78 65 34 98 12 23 54 34))
(emalgo N 0.2 20.0 10.0)

--=-=-=


-- 
Barry Fishman

--=-=-=--
From: Barry Fishman
Subject: Re: Lisp and Symbolic Integration
Date: 
Message-ID: <m3lk2afcew.fsf@barry_fishman.acm.org>
--=-=-=
Content-Type: text/plain; charset=utf-8
Content-Transfer-Encoding: 8bit


Although I made the changes to work with CLISP, I left the code
very non-portable.  I should really fix the still open issues.

Barry Fishman <·············@acm.org> writes:

> Francogrex <······@grex.org> writes:
>
>> On May 16, 12:16 am, Barry Fishman <·············@acm.org> wrote:
>>> So CLISP can actually handle larger numbers if they are expressed in
>>> long-float format.  If you want to keep you Lisp options open, you
>>> just need to avoid default conversions when needing large numbers.
>>
>> Thanks, now I see why this happens. Can we globally set all double to
>> long floats in clisp so that it behaves like GCL?
>
> You can do:
>
> (setf *read-default-float-format* 'long-float)

If you wish to compile your code for performance, one can not set just
set this value in an input file, since it will not be executed until
execution time, after your file has been read still interpreting
constant value strings as single precision.

You need to do this prior to compiling the file or include it within
an EVAL-WHEN so the compiler sees the change.

(eval-when (:compile-toplevel :load-toplevel :execute)
  (setf *read-default-float-format* 'long-float))

Also values should not be set prior to their being declared (either by a
LET or if they are to be global via DEFPARAMETER or DEFVAR.

So:
> ;;The mean function
> (defun mean(zazlist)
>   (setq sum (loop for x in zazlist sum x))
>   (setq m (/ sum (length zazlist))))

Really should be

(defun mean (zazlist)
  (let* ((sum (loop for x in zazlist sum x))
         (m (/ sum (length zazlist))))))

Or more simply:

(defun mean (zazlist)
  (/ (sum zazlist) (length zazlist)))

using your previously defined SUM function.

Also do:

;;Set eta (the tuning parameter)
(defparameter eta (list 1.0 1.0 1.0))

;;Set the counter.
(defparameter counter 0)

And later:

;; Run an example:
(defparameter N (list 23 24 78 65 34 98 12 23 54 34))

Rather than a SETQ (or SETF the more idiomatic assignment operator).
I changed COUNT to COUNTER because COUNT is part of the common lisp
package and is locked from being changed with more strict lisp systems.

In emalgo the varables OP1 OP2 OP3 and EZ are never declared and
so should be moved into a LET.

The update code is here:


--=-=-=
Content-Type: text/lisp
Content-Disposition: attachment; filename=tryme.lisp
Content-Description: EMALGO example code

;; The Expectation-Maximization (EM) algorithm

(eval-when (:compile-toplevel :load-toplevel :execute)
  (setf *read-default-float-format* 'long-float))

;; The while macro
(defmacro while (condition &rest body)
  (let ((var (gensym)))
    `(do ((,var nil (progn ,@body)))
	 ((null ,condition) ,var))))

;;The factorial function
(defun factorial (x) (if (zerop x) 1 (* x (factorial (- x 1)))))

;;The poisson PDF
(defun dpois (N L)
  (mapcar (lambda (x) (/ (* (expt L x) (exp (- L))) (factorial x)))
	  N))

;;The sum function
(defun sum (zazlist) (loop for x in zazlist sum x))

;;The mean function
(defun mean(zazlist)
  (/ (sum zazlist) (length zazlist)))

;;Set eta (the tuning parameter)
(defparameter eta (list 1.0 1.0 1.0))

;;Set the counter.
(defparameter counter 0)

;; The Expectation-Maximization (EM) algorithm function to find the
;; two means of the poisson mixtures (p2 & p3) and the mixing
;; probability p1.
(defun emalgo(N  p1 p2 p3)
  (while (or (> (nth 0 eta) 0.01)
	     (> (nth 1 eta) 0.01)
	     (> (nth 2 eta) 0.01))
    (incf counter)
    (let ((op1 p1) (op2 p2) (op3 p3))
      ;;The Expectation step
      (let ((EZ (mapcar #'/
			(mapcar (lambda (x) (* p1 x))
				(dpois N p2))
			(mapcar #'+
				(mapcar (lambda (x) (* p1 x))
					(dpois N p2))
				(mapcar (lambda (x) (* (- 1 p1) x))
				     (dpois N p3))))))

	;;The Maximization steps
	(setf p1 (mean EZ))
	(setf p2 (/ (sum (mapcar #'* EZ N )) (sum EZ)))
	(setf p3 (/ (sum (mapcar #'* (mapcar (lambda (x) (- 1.0 x)) EZ) N))
		    (sum (mapcar (lambda (x) (- 1.0 x))
				 EZ))))
	(setf eta (mapcar #'abs (mapcar #'- (list op1 op2 op3)
					(list p1 p2 p3))))

	;;Print out eta and the intermediate results of the p1,p2,p3 estimates.
	(print (list "eta:" eta))
	(print (list "count:" counter "p" (list p1 p2 p3)))))))

;; Run an example:
(defparameter N (list 23 24 78 65 34 98 12 23 54 34))
(emalgo N 0.2 20.0 10.0)

--=-=-=


It seem to work fine (without notes, warnings or errors) interpreted and
compiled for GCL, ECL, CLISP, SBCL, and CMUCL.  Those are the only lisp
systems I have available at the moment.

-- 
Barry Fishman

--=-=-=--
From: Pascal J. Bourguignon
Subject: Re: Lisp and Symbolic Integration
Date: 
Message-ID: <7c8wyaxl1f.fsf@pbourguignon.anevia.com>
Barry Fishman <·············@acm.org> writes:
> I have attached my own quick reformatted and changed code. (I made
> it into a unix style clisp script but you can just remove the initial #!
> line to LOAD it normally.)  It seems to give the same results as your
> program but with a few more digits displayed.

Note if you add:

(DEFUN EXECUTABLE-READER (A B C) (SYS::UNIX-EXECUTABLE-READER A B C))
(SET-DISPATCH-MACRO-CHARACTER #\# #\! (FUNCTION EXECUTABLE-READER))

to your ~/.clisprc, you'll be able to load these files without a hitch.

-- 
__Pascal Bourguignon__
From: Pascal J. Bourguignon
Subject: Re: Lisp and Symbolic Integration
Date: 
Message-ID: <7cd4nmzp4t.fsf@pbourguignon.anevia.com>
Francogrex <······@grex.org> writes:

> On May 15, 7:53�pm, Robert Dodier <·············@gmail.com> wrote:
>> On May 15, 9:32 am, Francogrex <······@grex.org> wrote:
>> > Maybe, but I don't think GCL is being sloppy, because R and S-plus
>> > work the same on this one as GCL and output exactly the same results.
>> > If GCL is sloppy then so is R (which seems strange).
>> That's beside the point; R is not an implementation of
>> Common Lisp, so it is not required to obey the same
>> rules about floating point numbers.
>
> Listen guys, you are the experts since you've been using lisp longer
> than I have. But compliant or not (though my version 2.6.7 is ansi
> compliant), GCL does what I need; when I type (exp 700) in GCL or
> Exp[700]//N in Mathematica I get what I need and I need large numbers
> and divisions by large numbers. Clisp isn't much help there for me.
> Just so that I know, is GCL looked down upon here?

Don't worry, that's the whole point of having various implementations
(even some that are not totally finickly compliant), that they may
have some differing behavior.

gcl and R behave the same because they just don't do anything.  They
rely on the underlying arithmetic processor.  You're just lucky that
you're running both programs on the same processor, that implements
the same IEEE754 floating point arithmetics.  


Anyways, for your purpose, gcl might be a good choice, given that
maxima is well ported on gcl.  It also works well on clisp, sbcl, etc,
so you can try it with all these implementations, and commit to one
later.

http://maxima.sourceforge.net/lisp.html

-- 
__Pascal Bourguignon__