[scilab-Users] advice on least square method

[ray joseph]

Paul,

If you have a look at http://www.mathworks.com/moler/leastsquares.pdf you will find not only a vectorized method, but considerations on what to do when the data/experient is ill conditioned.  Section 5.5 shows the easy way to calculate the matrix of sums of powers.  One consideration is to not calculate the determinant to see if there is a solution.  This calc is very expensive.  Let the algorythm detect singularity.

Ray
  ----- Original Message ----- 
  From: Carrico, Paul 
  To: users at lists.scilab.org 
  Sent: Friday, August 27, 2010 5:04 AM
  Subject: [scilab-Users] advice on least square method


  Dear all,



  The current email is intend as (an) advice(s) request regarding the least square method calculations.



  I'm currently using the function herebellow to calculate the coefficients of a 2nd order polynomial using thousands of experimental data (furthermore the function is inserted in a loop => thousands of loop = billiards of calculations).



  So the main time cost is these calculation ...





  I'm working to improve the code to be faster and in parallele I would like to know if these a more relevant instruction in Scilab that performe the calculation faster (maybe more accurate) ?



  I had a look in the leastsq instruction ... but the later is strongly influenced by the initial parameters guess X0 (isn't it).





  Any advice is welcome



  Regards / have a good WE



  Paul



  PS : A give a basic exampel in attached

  PS 2 : additional question (but i'm a newby : calculations to be parallelized ???)







  ###################################################"





  ##########################################################

  function R_square = fct_MMC(p,experience)

   

  // Nota :

  // p = polynomial degree

  // experience = experimental data





  // step 1 : "linear" least square method

  // matrix size
  [n,c] = size(experience);

   

  // calculation of  S
  S = zeros(2*p,1);
  for k = 1 : (2*p)
   for i = 1 : n
    S(k) = S(k) + experience(i,1)^k;
   end
  end

   

  // calculation of W
  W = zeros(p+1,1);
  for k = 1 : (p+1)
   for i = 1 : n
    W(k) = W(k) + experience(i,2)*experience(i,1)^(k-1);
   end
  end

   

  // calculation of A = S matrix
  A = zeros(p+1,p+1);
  for i = 1 : (p+1)
   for j = 1 : (p+1)
    if (i == 1 & j == 1) then
     A(1,1) = n;
    else
     A(i,j) = S(i+j-2);
    end
   end
  end

   


  // calculation of  B = W matrix
  B = zeros(p+1,1);
  for j = 1 : (p+1)
   B(j) = W(j);
  end

   


  // calculation of the coefficients
  if (det(A) == 0) then
   printf(' #######################################################################################\n')
   printf(' ####                                              Singular matrix => no solution                                                 ####\n')
   printf(' #######################################################################################\n')

   

  else
   //C = inv(A)*B;

  C = lsq(A,B);   // a bit faster than the previous calculation / same results
   // calculation of R_square
   // formula : R_square = 1 - ( sum(F_fit - F_measured)^2 / sum(F_mean - F_measured )
   F_fit = zeros(n,1);
   R_square = 0;
   F_mean = 0;
   
   for i = 1 : n
    // Calcul de F_fit
     for k = 1 : (p+1)
      F_fit(i) = F_fit(i) + C(k)*experience(i,1)^(k-1) ;
     end
   end

   

   calculation of  F_mean
   F_mean = mean(experience(:,2));

   

   // SSE = sum (F_fit - F_mean)
   // SSTE = sum (F_mean - F_measured)
   SSE = 0;
   SST = 0;



   for i = 1 : n
    SSE = SSE + (F_fit(i) - experience(i,2))^2 ;
    SST = SST + (F_mean - experience(i,2))^2 ;
   end

   

   R_square = 1 - (SSE / SST);

   

   

  end

   

  clear A
  clear B
  clear S
  clear W
  clear SSE
  clear SST

   


  endfunction

   

   

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


Le présent mail et ses pièces jointes sont confidentiels et destinés à la personne ou aux personnes visée(s) ci-dessus. Si vous avez reçu cet e-mail par erreur, veuillez contacter immédiatement l'expéditeur et effacer le message de votre système. Toute divulgation, copie ou distribution de cet e-mail est strictement interdite.

This email and any files transmitted with it are confidential and intended solely for the use of the individual or entity to whom they are addressed. If you have received this email in error, please contact the sender and delete the email from your system. If you are not the named addressee you should not disseminate, distribute or copy this email.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20100827/c1d2911f/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: matrix.JPG
Type: image/jpeg
Size: 24546 bytes
Desc: not available
URL: <https://lists.scilab.org/pipermail/users/attachments/20100827/c1d2911f/attachment.jpe>