Re: FW: lspoly didn't wok when the number of data points equal to the number of coefficients

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Fri Jan 03 2014 - 13:25:54 MST

Attached is the SLATEC fortran code and a sample test script
that calls the fortran subroutine. You will have to change the
NCL driver script.

%> WRAPIT slatec_DPOLFT.f

Please note: This is a double precision version of the
fortran code. The data passed to it must be dpoble
precision. Hence, if x and y are type float

  xd = xf*1d0 ; new array in double precisiom
  yd = yf*1d0

LSPOLY::POLFT (M, xd, yd, wgt ...)

===

Please read the Mini-Language Manual for calling fortran code

http://www.ncl.ucar.edu/Document/Manuals/

Good luck

On 01/02/2014 07:17 AM, Richard Liang wrote:
> Hi, Dennis
>
> I could not wait the next version of NCL since my project is urgent.
> Can I use the lspoly_ncl.f you sent to me in NCL? How to use it ?
> Do you have the slatec_lspoly function in NCL
>
> Thanks
>
>
>
> > Date: Wed, 1 Jan 2014 09:13:04 -0700
> > From: shea@ucar.edu
> > To: richardljy@hotmail.co.uk; ncl-talk@ucar.edu
> > Subject: Re: FW: lspoly didn't wok when the number of data
> points equal to the number of coefficients
> >
> > The lspoly used thru NCL v6.1.2 has some known limitations.
> > The next release (v6.2.0) will hopefully contain an updated version of
> > lspoly. Here called 'slatec_lspoly'. The following test script
> > and the results illustrates the issues.
> >
> > A comparison to R results follows.
> >
> > ==> NCL script
> >
> > x = (/-4.5, -3.2, -1.4/)
> > y = (/ 0.7, 2.3, 3.8/)
> >
> > n = dimsizes(y)
> >
> > print(" ")
> > print("===========> 6.1.2 <===========")
> > print(" ")
> >
> > n1 = n-1 ; max # of coef returned by lspoly
> > ; n1=2
> > c = lspoly(x,y, 1, n1) ; up thru 6.1.2
> > cc = slatec_lspoly(x,y, 1, n1) ; 6.2.0
> >
> > print(c)
> > print(cc)
> >
> > cx = lspoly(x,y, 1, n) ; n=3 will fail
> > print(cx)
> >
> > print(" ")
> > print("===========> 6.2.0 <===========")
> > print(" ")
> >
> > cc3 = slatec_lspoly(x,y, 1, n) ; n=3 here
> > print(cc3)
> >
> > *************************************************
> >
> > (0)
> > (0) ===========> 6.1.2 <===========
> > (0)
> >
> > Variable: c
> > Type: float
> > Total Size: 8 bytes
> > 2 values
> > Number of Dimensions: 1
> > Dimensions and sizes: [2]
> > Coordinates:
> > (0) 5.268707
> > (1) 0.9896836
> >
> > Variable: cc
> > Type: float
> > Total Size: 8 bytes
> > 2 values
> > Number of Dimensions: 1
> > Dimensions and sizes: [2]
> > Coordinates:
> > (0) 5.268707
> > (1) 0.9896836
> >
> > Variable: cx <=== fail
> > Type: float
> > Total Size: 12 bytes
> > 3 values
> > Number of Dimensions: 1
> > Dimensions and sizes: [3]
> > Coordinates:
> > (0) 1e+20
> > (1) 1e+20
> > (2) 1e+20
> >
> > (0)
> > (0) ===========> 6.2.0 <===========
> > (0)
> >
> > Variable: cc3 match R results
> > Type: float
> > Total Size: 12 bytes
> > 3 values
> > Number of Dimensions: 1
> > Dimensions and sizes: [3]
> > Coordinates:
> > (0) 4.392307
> > (1) 0.2435896
> > (2) -0.1282052
> >
> >
> >
> > ***********************************************
> > ------ Using R --------------------------------
> > ************************************************
> >
> > > x <- c(-4.5, -3.2, -1.4)
> > > y <- c(0.7, 2.3, 3.8)
> > > n <- 1 # degree of polynomial
> > > result <- lm( y ~ poly( x, n, raw=T))
> >
> > > result
> >
> > Call:
> > lm(formula = y ~ poly(x, n, raw = T))
> >
> > Coefficients:
> > (Intercept) poly(x, n, raw = T)
> > 5.2687 0.9897
> >
> > ---------------------
> >
> > > n <- 2 # degree of polynomial
> > > result <- lm( y ~ poly(x,n, raw=T))
> > > result
> >
> > Call:
> > lm(formula = y ~ poly(x, n, raw = T))
> >
> > Coefficients:
> > (Intercept) poly(x, n, raw = T)1 poly(x, n, raw = T)2
> > 4.3923 0.2436 -0.1282
> >
> > On 12/31/13, 8:18 PM, Richard Liang wrote:
> > > Thanks. Dennis. But the problem remain: Mathematically, if I have 3
> points, I can make an second order polynomial to fit these points , in
> which I should have 3 coefficients.
> > > Thanks
> > >
> > >> Date: Mon, 30 Dec 2013 09:35:37 -0700
> > >> From: shea@ucar.edu
> > >> To: richardljy@hotmail.co.uk; ncl-talk@ucar.edu
> > >> Subject: Re: FW: lspoly didn't wok when the number of
> data points equal to the number of coefficients
> > >>
> > >> The Description section documentation was not correct. It has been
> changed.
> > >>
> > >> On 12/29/13, 9:53 AM, Richard Liang wrote:
> > >>> Thanks, Dennis
> > >>> The document said "It is necessary that the number of data
> points) be greater than or equal to n(the number of coefficients)."
> > >>> As I understand , for example, if I have 3 points, the max number
> of coefficients can be 3.
> > >>> Mathematically, if I have 3 points, I can make an secend order
> polynomial to fit these points , in which I should have 3 coefficients.
> > >>>
> > >>> Best
> > >>>
> > >>>> Date: Sat, 28 Dec 2013 08:26:49 -0700
> > >>>> From: shea@ucar.edu
> > >>>> To: richardljy@hotmail.co.uk; ncl-talk@ucar.edu
> > >>>> Subject: Re: FW: lspoly didn't wok when the number of
> data points equal to the number of coefficients
> > >>>>
> > >>>> As noted in the documentation:
> > >>>>
> > >>>> n
> > >>>>
> > >>>> The number of coefficients desired (i.e., n-1 will be the degree
> of the
> > >>>> polynomial). Due to the method used, n should be less than or
> equal to
> > >>>> five.
> > >>>>
> > >>>> ===
> > >>>>
> > >>>> In your case with 3 points, the max number of coefficients would
> be 2
> > >>>>
> > >>>>
> > >>>> On 12/28/13, 5:01 AM, Richard Liang wrote:
> > >>>>>
> > >>>>>
> > >>>>> From: richardljy@hotmail.co.uk
> > >>>>> To: ncl-talk@ucar.edu
> > >>>>> Subject: lspoly didn't wok when the number of data points equal
> to the number of coefficients
> > >>>>> Date: Sat, 21 Dec 2013 15:45:08 -0500
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>> lspoly didn't wok when the number of data points equal to the
> number of coefficientsI did the test
> > >>>>> x = (/-4.5, -3.2, -1.4/)y = (/ 0.7, 2.3, 3.8/)n = 3 c =
> lspoly(x,y, 1, n) print(c)c will be(0) 1e+20(1) 1e+20(2) 1e+20
> > >>>>>
> > >>>>>
> > >>>>> But if n=2 c will be(0) 5.268707(1) 0.9896836
> > >>>>> I also did the test in the 4 element array, when n=4, it didn't
> work. only when n is less than the number of data points
> > >>>>> Any body help me? thanks
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>>
> > >>>>> _______________________________________________
> > >>>>> ncl-talk mailing list
> > >>>>> List instructions, subscriber options, unsubscribe:
> > >>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> > >>>>>
> > >>>
> > >>>
> > >
> > >

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk

Received on Fri Jan 03 13:26:06 2014

This archive was generated by hypermail 2.1.8 : Mon Jan 06 2014 - 13:09:54 MST