BookmarkSubscribeSubscribe to RSS Feed

Re: Problem extracting derivative of Spline fit in JSL

yvesprairie

Community Trekker

Joined:

Nov 25, 2014

HI everyone,

 

I am having a curious problem I cannot figure out. I have 2 variables (say Y and X) for which I want to extract the derivative of a smooth function linking the two. I can do it "manually" if my two variable are in a table. I simply fit a smooth spline, save the predicted formula, go in the formula editor and use the derivative option in the editor itself. Now I want to do the same but in an automated fashion because I will have dozens of such fitting to do so have reverted to using JSL. As far as I can tell, the only to get there in JSL is first to extract the spline coefficients as

 

splcoef=Spline Coef(X, Y, 1);

 

In the particular example that I am using this matrix would be:

 

[0.5 0.997009092716414 0.00052834644305232 0 0.0000661345992050551,

1 0.997281532762841 0.000577947392456002 0.0000992018988076818 0.0000439565382573196,

1.5 0.997600801501053 0.000710116694956763 0.000165136706193727

-0.0000523992394403174,

2 0.99799059412015 0.000835953971570102 0.0000865378470331725 -0.0000897493429406088,

2.5 0.998418986899826 0.00085517981139797 -0.0000480861673778753

-0.0000494911432825247,

3 0.99882836887077 0.000769975286558009 -0.000122322882301737 -0.0000204292892022263,

3.5 0.999180222132323 0.000632330437354731 -0.000152966816105107

-0.0000100612562586075,

4 0.999456887989942 0.000471817679055371 -0.000168058700493033 0.0000251846412034386,

4.5 0.999653930234497 0.000322647459465303 -0.000130281738687837

0.0000199230550287456,

5 0.999785173911436 0.000207308012048703 -0.000100397156144689 0.0000256287090160626,

5.5 0.999866932217051 0.000126132387666307 -0.0000619540926205567

0.0000114058197541642,

6 0.999915935615199 0.00007273265986134 -0.0000448453629892932 0.0000100640302564247,

7 0.999953886942327 0.0000132340246520343 -0.0000146532722199889

0.0000064471411305935,

7.5 0.999957646529239 0.0000034161082800262 -0.000004982560524089

0.0000026947675920361,

8 0.999958445789198 0.0000004546234498223 -0.0000009404091360308

0.0000012459427765851,

8.5 0.999958593741486 0.000000448671396435 0.0000009285050288488

-0.0000002521994682104,

9 0.999959018678507 0.0000011880268241599 0.0000005502058265328

-0.0000001515710315195,

9.5 0.999959731296997 0.0000016245543768191 0.0000003228492792533

-0.0000002152328528353,

10 0.999960597382399 0 0 0]

 

I can find the predicted values from the spline as:

 

predvalues=Spline Eval(X,splcoef);

 

to obtain the derivative, I tried to use:

 

dervalues=eval(derivative(spline eval(X, splcoef),X));

 

but this gives a unique value of zero (0).

 

However, if I replace the matrix reference splcoef by its actual contents in the spline eval statement then it works fine!!

 

dervalues=eval(derivative(spline eval(X,

[0.5 0.997009092716414 0.00052834644305232 0 0.0000661345992050551,

1 0.997281532762841 0.000577947392456002 0.0000992018988076818 0.0000439565382573196,

1.5 0.997600801501053 0.000710116694956763 0.000165136706193727

-0.0000523992394403174,

2 0.99799059412015 0.000835953971570102 0.0000865378470331725 -0.0000897493429406088,

2.5 0.998418986899826 0.00085517981139797 -0.0000480861673778753

-0.0000494911432825247,

3 0.99882836887077 0.000769975286558009 -0.000122322882301737 -0.0000204292892022263,

3.5 0.999180222132323 0.000632330437354731 -0.000152966816105107

-0.0000100612562586075,

4 0.999456887989942 0.000471817679055371 -0.000168058700493033 0.0000251846412034386,

4.5 0.999653930234497 0.000322647459465303 -0.000130281738687837

0.0000199230550287456,

5 0.999785173911436 0.000207308012048703 -0.000100397156144689 0.0000256287090160626,

5.5 0.999866932217051 0.000126132387666307 -0.0000619540926205567

0.0000114058197541642,

6 0.999915935615199 0.00007273265986134 -0.0000448453629892932 0.0000100640302564247,

7 0.999953886942327 0.0000132340246520343 -0.0000146532722199889

0.0000064471411305935,

7.5 0.999957646529239 0.0000034161082800262 -0.000004982560524089

0.0000026947675920361,

8 0.999958445789198 0.0000004546234498223 -0.0000009404091360308

0.0000012459427765851,

8.5 0.999958593741486 0.000000448671396435 0.0000009285050288488

-0.0000002521994682104,

9 0.999959018678507 0.0000011880268241599 0.0000005502058265328

-0.0000001515710315195,

9.5 0.999959731296997 0.0000016245543768191 0.0000003228492792533

-0.0000002152328528353,

10 0.999960597382399 0 0 0]

),X));

 

All I have done is replace the matrix name by its contents. This does not seem like normal behavior.

 

Any ideas how I get it to work??

 

Many thanks, Yves P.

 

 

 

 

1 ACCEPTED SOLUTION

Accepted Solutions
Highlighted
ms

Super User

Joined:

Jun 23, 2011

Solution

Try the Eval( Evalexpr( ) ) cunstruct to force preevaluation of the matrix.

 

Names Default To Here(1);
X = 1 :: 100;
Y = Normal Distribution(X, 50, 25);
splcoef = Spline Coef(X, Y, 1);
predvalues = Spline Eval(X, Spline Coef(X, Y, 1));
dervalues = Eval(Eval Expr(Derivative(Spline Eval(X, Expr(splcoef)), X)));

// Plot to see if it works
dt = New Table("test", add rows(100));
dt << New Column("X", values(X)) << New Column("Y", values(predvalues)) <<
New Column("dY", values(dervalues));

dt << Graph Builder(
    Variables(X(:X), Y(:Y), Y(:dY)),
    Elements(Position(1, 1), Points(X, Y)),
    Elements(Position(1, 2), Points(X, Y))
);
2 REPLIES
Highlighted
ms

Super User

Joined:

Jun 23, 2011

Solution

Try the Eval( Evalexpr( ) ) cunstruct to force preevaluation of the matrix.

 

Names Default To Here(1);
X = 1 :: 100;
Y = Normal Distribution(X, 50, 25);
splcoef = Spline Coef(X, Y, 1);
predvalues = Spline Eval(X, Spline Coef(X, Y, 1));
dervalues = Eval(Eval Expr(Derivative(Spline Eval(X, Expr(splcoef)), X)));

// Plot to see if it works
dt = New Table("test", add rows(100));
dt << New Column("X", values(X)) << New Column("Y", values(predvalues)) <<
New Column("dY", values(dervalues));

dt << Graph Builder(
    Variables(X(:X), Y(:Y), Y(:dY)),
    Elements(Position(1, 1), Points(X, Y)),
    Elements(Position(1, 2), Points(X, Y))
);
yvesprairie

Community Trekker

Joined:

Nov 25, 2014

Thank you so much, you just saved me hours of work!!! Sincerely, Yves