Topic: [Challenge] least square polynomial

[dgpuertas]

The challenge is to find the polynomial of 'n' degree which leaves the minimum distance to a given points.
I use the following method:


And this is my code:

Code - Auto/Visual Lisp: [Select]

1. (defun coeficientes_n_grado_minimos_cuadrados (lis pot / sublist n x y lsx lsy mat)
2.    
3.   (setq n (length lis)
4.         x (mapcar (function car) lis)
5.         y (mapcar (function cadr) lis)
6.         lsx (cons n
7.                   (mapcar (function (lambda (ex)
8.                         (apply (function +)
9.                                (mapcar (function (lambda (a) (expt a ex))) x))))
10.                           (iseq 1 (* pot 2))))
11.         lsy (cons (apply (function +) y)
12.                   (mapcar (function (lambda (ex)
13.                         (apply (function +) (mapcar (function (lambda (a b) (* (expt a ex) b))) x y))))
14.                           (iseq 1 pot)))
15.         mat (mapcar (function (lambda (n)
16.                 (mapcar (function (lambda (a) (nth a lsx))) (iseq n (+ n pot)))))
17.                     (iseq 0 pot))
18.   )
19.  
20. (if (> n pot)
21.     (mxv (inverse_matrix mat)
22.        lsy))
23. )
24.  
25. (defun iseq (i j / ls)
26.   (repeat (1+ (- j i))
27.     (setq ls (cons j ls)
28.           j (1- j)))
29. ls)
30.  
31.  

I use inverse_matrix and dot vector matrix function to resove the ecuations.

The function 'coeficientes_n_grado_minimos_cuadrados' needs a lists of 2Dpoints and a degree of the polynomial and returns the coeficients of the polynomial.

Example:

Code - Auto/Visual Lisp: [Select]

1. (coeficientes_n_grado_minimos_cuadrados
2. '('(73.5017 32.4285 0.0) ('110.745 0.519587 0.0) ('172.952 0.519587 0.0) ('204.058 23.4271 0.0) ('224.933 60.2435 0.0) ('285.916 69.6531 0.0))
3. 3)
4.  
5.  
6. (300.845 -5.69764 0.0322512 -5.29193e-005)
7.  

represent polynomial 3 degree:
300.845  -5.69764 x + 0.0322512x² -5.29193e-005x³


Thanks and sorry about my english

[chlh_jd]

Nice codes , Dgpuertas  .
Thanks you to share it . 

[dgpuertas]

Thanks a lot chlh_jd

The code works with any degree, since a straight line (with 1degree) until n degree (with n less than number of points)


Thanks and sorry about my english

[chlh_jd]

Thanks a lot chlh_jd
The code works with any degree, since a straight line (with 1degree) until n degree (with n less than number of points)

You're welcome !

I think , If the Pointset faraway the origin of coordinates it will take get an undesirable result , So we must translate the given pointset into a Assumed coordinate system which X-axis does not contain any point of the given pointset .
I use your function like following codes .

Code: [Select]

(defun c:Test  (/ ptlst pt Xlst Ylst i k n lst x0 xn div mid ss s1
    sslst midlst ent1 ent2 pt0 a b str osm)
  ;; Use dgpuertas's least-squares-solution function to Regress Curve
  ;; by GSLS(SS) 2013-2-21
  (setq pt
(getpoint
   "\nGet the first regress points (RightButton into Select all regress points)>>:"))
  (while pt
    (setq ptlst (cons pt ptlst)
  pt nil
  pt (getpoint "\nContinue selecting the next point>>:")))
  (if (not ptlst)
    (progn
      (prompt "\nSelect all regress points or cross-lines:")
      (setvar "nomutt" 1)
      (setq ss (ssget '((0 . "LINE,POINT"))))
      (setvar "nomutt" 0)
      (if ss
(progn
  (setq i     0
sslst nil)
  (while (setq s1 (ssname ss i))
    (setq sslst (cons s1 sslst))
    (setq i (1+ i)))
  (foreach a  sslst
    (foreach b sslst
      (setq ent1 (entget a)
    ent2 (entget b))
      (if (= (cdr (assoc 0 ent1)) "POINT")
(setq ptlst (cons (cdr (assoc 10 ent1)) ptlst))
(if (= (cdr (assoc 0 ent2)) "POINT")
  (setq ptlst
(cons (cdr (assoc 10 ent2)) ptlst))
  (if (setq pt (inters (cdr (assoc 10 ent1))
       (cdr (assoc 11 ent1))
       (cdr (assoc 10 ent2))
       (cdr (assoc 11 ent2))))
    (setq ptlst (cons pt ptlst)))))))
  (setq midlst nil)
  (foreach a  ptlst
    (if (not (member a midlst))
      (setq midlst (cons a midlst))))
  (setq ptlst midlst))
(princ "\nNo selection !")
)
      )
    )
  (if ptlst
    (progn
      (setq ptlst (vl-sort
    ptlst
    (function (lambda (e1 e2) (< (car e1) (car e2)))))
    pt0   (midpts ptlst)
    ptlst
  (mapcar (function
    (lambda (x)
      (gsls-XY->AB x pt0 0)))
  ptlst))
      (or (setq k
(getint
   (strcat
     "\nType in RegressEquation Degree (Should be not more than the half of the pointset "
     (rtos (/ (length ptlst) 2) 2 0)
     ", and must be less than the total of the pointset <"
     (rtos (length ptlst) 2 0)
     ") <2>:")))
  (setq k 2))
      (if (setq lst (2dcurve:least-squares-solution ptlst k))
(progn
  (setq str " "
i   0)
  (repeat (1+ k)
    (setq str
   (strcat str
   (strcat (if (minusp (nth i lst))
     ""
     "+")
   (rtos (nth i lst) 1 6)
   (if (= i 0)
     ""
     (strcat "X^" (rtos i 2 0))))))
    (setq i (1+ i)))
  (setq str (strcat "\nRelatively pt0 "
    (vl-princ-to-string pt0)
    " curve equation is : Y="
    (vl-string-trim "+" str)))
  (princ str)
  (or (setq n
     (getint
       "\nType in the devided number of the X-coordinate segment to draw the regress curve <100>:"))
      (setq n 100))
  (setq x0   (caar ptlst)
xn   (car (last ptlst))
div  (/ (- xn x0) n)
i    0
Xlst nil)
  (while (<= i n)
    (setq Xlst (cons (+ x0 (* i div)) Xlst)
  i    (1+ i)))
  (setq Ylst nil)
  (foreach a  Xlst
    (setq i   0
  mid 0)
    (repeat (1+ k)
      (setq mid (+ mid
   (* (if (= a 0)
0.0
(expt a i))
      (nth i lst)))
    i (1+ i)))
    (setq Ylst (cons mid Ylst)))
  (setq ptlst nil
i     0)
  (while (setq Xi (nth i Xlst))
    (setq Yi (nth (- n i) Ylst)
  ptlst (cons (list Xi Yi) ptlst)
  i (1+ i)))
  (setq ptlst
(mapcar (function
   (lambda (x)
     (gsls-AB->XY x pt0 0)))
ptlst))
  (setq osm (getvar "OSMODE"))
  (setvar "OSMODE" 0)
  (command "_.pline")
  (apply (function command) ptlst)
  (command "")
  (setvar "OSMODE" osm))
(princ
  "\nError:regression solving error , please check pointset and parameters ")))
    (princ "\nNo Selected Pointset !")
    )
  (princ "\nNonlinear regression program completed.")
  (princ)
  )