Carl Love

Carl Love

28100 Reputation

25 Badges

13 years, 105 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@Markiyan Hirnyk 

When I said that my fit had lower residuals, I was using the fit that treated the first four points as outliers (my Sol_Carl5), and I compared it to doing the equivalent with DirectSearch:-DataFit, using all nine fitmethods and taking the minimal residuals. Done this way, my method has less than half the residuals.

But, if I use all four methods crossed with all nine fitmethods (for a total of 36 calls to DataFit), I can make a small improvement over my method. Here's the worksheet:

 

restart:

TW:= <<2.00E-05|0.00723>,<4.86E-05|0.033171>,<7.31E-05|0.047349>,<9.77E-05|0.057941>,<0.000122|0.066835>,<0.000147|0.074711>,<0.000171|0.081864>,<0.000196|0.088463>,<0.000221|0.094621>,<0.000245|0.10042>,<0.00027|0.10591>,<0.000295|0.11114>,<0.000319|0.11614>,<0.000344|0.12091>,<0.000368|0.12544>,<0.000393|0.12977>,<0.000418|0.13391>,<0.000442|0.13788>,<0.000467|0.1417>,<0.000492|0.14539>,<0.000516|0.14896>,<0.000541|0.15244>,<0.000565|0.15582>,<0.00059|0.15911>,<0.000615|0.16231>,<0.000639|0.16544>,<0.000664|0.1685>,<0.000689|0.17148>,<0.000713|0.1744>,<0.000738|0.17725>,<0.000762|0.18005>,<0.000787|0.18279>,<0.000812|0.18548>,<0.000836|0.18811>,<0.000861|0.1907>,<0.000885|0.19324>,<0.00091|0.19574>,<0.000935|0.19819>,<0.000959|0.2006>,<0.000984|0.20297>,<0.001008|0.20531>,<0.001033|0.20761>,<0.001058|0.20987>,<0.001082|0.2121>,<0.001107|0.2143>,<0.001131|0.21647>,<0.001156|0.2186>,<0.00118|0.22071>,<0.001205|0.22279>,<0.00123|0.22484>,<0.001254|0.22687>,<0.001279|0.22886>,<0.001303|0.23083>,<0.001328|0.23278>,<0.001352|0.2347>,<0.001377|0.2366>,<0.001401|0.23847>,<0.001425|0.24032>,<0.00145|0.24215>,<0.001474|0.24396>,<0.001498|0.24575>,<0.001523|0.24751>,<0.001547|0.24926>,<0.001571|0.25099>,<0.001596|0.25277>,<0.00162|0.25433>,<0.001644|0.25571>,<0.001668|0.25708>,<0.001693|0.25843>,<0.001717|0.25978>,<0.001741|0.26112>,<0.001765|0.26245>,<0.001789|0.26377>,<0.001813|0.26508>,<0.001838|0.26638>,<0.001862|0.26767>,<0.001886|0.26895>,<0.00191|0.27022>,<0.001934|0.27148>,<0.001958|0.27273>,<0.001982|0.27396>,<0.002006|0.27519>,<0.00203|0.27641>,<0.002054|0.27762>,<0.002078|0.27881>,<0.002102|0.28>,<0.002126|0.28118>,<0.00215|0.28235>,<0.002174|0.28352>,<0.002198|0.28467>,<0.002222|0.28582>,<0.002246|0.28695>,<0.00227|0.28808>,<0.002294|0.2892>,<0.002318|0.29032>,<0.002342|0.29142>,<0.002366|0.29252>,<0.00239|0.29361>,<0.002414|0.29469>,<0.002438|0.29577>,<0.002462|0.29684>,<0.002486|0.2979>,<0.002509|0.29895>,<0.002533|0.3>,<0.002557|0.30104>,<0.002581|0.30207>,<0.002605|0.3031>,<0.002629|0.30412>,<0.002653|0.30513>,<0.002677|0.30614>,<0.002701|0.30714>,<0.002724|0.30813>,<0.002748|0.30912>,<0.002772|0.3101>,<0.002796|0.31107>,<0.00282|0.31204>,<0.002844|0.313>,<0.002868|0.31396>,<0.002891|0.31491>,<0.002915|0.31586>>:

RecenteredPowerFit:= proc(TW::Matrix, b_range::range(realcons))

local

     T:= TW[..,1], W:= TW[..,2],

     b:= Optimization:-NLPSolve(

          b-> Statistics:-PowerFit(T-~b, W, output= residualsumofsquares),

          b_range, method= branchandbound, objectivetarget= 0

     )[2][1],

     Res:= Statistics:-PowerFit(T-~b, W)

;

     ['A' = Res[1], ':-b' = b, 'c'= Res[2]]

end proc:

 

Sol_Carl5:= RecenteredPowerFit(TW[5.., ..], -1..TW[5,1] - 1e-7);

[A = HFloat(3.682236898705571), b = HFloat(6.28308061404132e-5), c = HFloat(0.4158452235037357)]

model:= A*(t-b)^c:

FitMethods:= [lms, lad, lts, lws, minimax, median, quantile, trimean, lfad]:

T5:= TW[5.., 1]:  W5:= TW[5.., 2]:
Min:= infinity:
for f in FitMethods do
     Sol_DS[f]:= DirectSearch:-DataFit(
          model, T5, W5, t,
          fitmethod= f, evaluationlimit= 30000, usewarning= false
     );
     V:= add(abs((eval(eval(model, Sol_DS[f][2]), t= T5[k]) - W5[k])^2), k= 1..op(1,T5));
     if V < Min then  (Min,Min_f):= V,f  end if  
end do:

RSS_DS:= Min;

HFloat(0.005773814784137672)

Min_f;

trimean

RSS_Carl:= add((eval(eval(model, Sol_Carl5), t= T5[k]) - W5[k])^2, k= 1..op(1,T5));

HFloat(8.687836999527558e-4)

Now try all methods and all fitmethods of  DataFit.

Methods:= [cdos, brent, powell, quadratic]:

Min:= infinity:
for f in FitMethods do  for m in Methods do
     Sol_DS[f,m]:= DirectSearch:-DataFit(
          model, T5, W5, t,
          fitmethod= f, evaluationlimit= 30000, usewarning= false,
          method= m, strategy= globalsearch
     );
     V:= add(abs((eval(eval(model, Sol_DS[f,m][2]), t= T5[k]) - W5[k])^2), k= 1..op(1,T5));
     if V < Min then  (Min,Min_f,Min_m):= V,f,m  end if  
end do  end do:

 

RSS_DS:= Min;

HFloat(6.044545213233039e-4)

Min_f, Min_m;

lms, quadratic

 

 

Download RecenteredPowerFit_vs_DirectSeach.mw

@Markiyan Hirnyk I'm not trying to blame anybody! I was just trying to tell you that both issues were caused by the same mistake. Geez, I can't even say that you're right without you finding some way to criticize!

@capea 

You didn't execute the expandoff command! That is nescessary to prevent simplification of residues. Also, don't forget to restart before the whole process.

 

@Markiyan Hirnyk You're absolutely right. I made a mistake. I will correct the Answer now.

Regarding the expand portion, you're experiencing the same mistake. After I make the correction to the Answer, try expand(Psi(3)) before and after the expandoff commands.

@MDD Try using the Maple help: Enter ?Psi at the Maple prompt.

Psi(z) = diff(ln(GAMMA(z)), z)

@MDD You simply cut and paste and execute that code that Roman gave, then 

with(PolynomialIdeals):
EliminationIdeal(<x>,{a,b,c});

@Markiyan Hirnyk 

By what criterion are you saying that the result produced by DirectSearch:-DataFit is better than that produced by my RecenteredPowerFit (which calls Statistics:-PowerFit)? You can't directly compare the residuals as reported by the various fitmethods because they are each computed by their own formula. If you compare the sum of the squares of the residuals, my fit is several times better than yours. It is also clear from the plots that my fit is better than yours. (My fit is the red curve in my plot; ignore the blue curve.)

Any result with > 1 is completely ridiculous, not just "not so good". Try plotting them. Of the nine fitmethods, seven produce such results for me (on the first run).

Furthermore, the results of DirectSearch:-DataFit are not entirely reproducible. Rerunning the exact same call (of course, without restart) often gives a completely different result, often changing a decent result into a completely ridiculuous one.

The vote up may be from the OP and may be due to the fact that I worked back and forth with her over three days and that I produced code that she was able to execute and with which she was able to get what she considered to be "perfect" results.

 

@Jerome BENOIT I also don't care what the shown material is called. I just didn't want other readers to get the impression that it was ever possible to see the complete code of a module---not by eval or by any other means. Your usage of "body" was likely to give that false impression.

Here's a suggestion: stabilize the plot to the unit circle by dividing by v; get rid of the view option. This makes animations that vary n or m much easier. So here's my version, a modification of your cycler from Reply "without t", not the one with the circles. I also added _rest to the plot command to make it easier to pass in plot options.

cycler := proc(k, p, m, n, T) local expr, t, u, up, v;
     u := exp(k*I*t);
     up := exp(-k*I*t*p);
     expr := exp(I*t) * (1 - u/m + I*up/n);
     v := 1 + abs(1/n) + abs(1/m);
     plots:-complexplot( expr/v, t = 0 .. T, axes = none, _rest);
end proc:

And here's an animation varying n. Note that I vary it exponentially.

plots:-animate(
     cycler, [5, 3, 2, 10^t, 2*Pi, color= purple], t= -2..1.5,
     frames= 100, paraminfo= false
);

And here's one varying m, also exponentially:

plots:-animate(
     cycler, [5, 3, 10^t, 3, 2*Pi, color= orange], t= -2..2.5,
     frames= 100, paraminfo= false
);

@Jerome BENOIT 

No version of Maple has ever shown a module's body, which is the executable code after the declarations. Indeed, that code is discarded after it's executed; it's not stored in the module. What eval shows (or used to show) is the module's header or shell, which includes the description, the options, and the names (but not the values, if any) of the locals and exports.

@Alejandro Jakubi What if you set interface(verboseproc= 2)?

@MDD 

module is a Maple data structure similar to a procedure that allows for a special type of local variable called an export which can be accessed outside the module. It is the primary structure in Maple for object-oriented programming, similar to a class in C++ or Java. Often, the exports of a module are procedures which allow controlled access to the locals. See ?module. Sometimes a module is used primarily or even exclusively as a container object like a table. See ?record.

This has nothing to do with the algebraic structure known as a module.

@Haley_VSRC 

Markiyan is using the DirectSearch package, which must be downloaded from the Maple Applications Center. It's a very good package with many different algorithms for nonlinear fitting. But, as you can see, several of those algorithms return "crazy" results with > 1.

@Haley_VSRC 

Here's how to use my procedure. In the final graph, note that there's a trend in the rightmost quarter or so of the data that is not explained by either curve.


restart:

TW:= <<2.00E-05|0.00723>,<4.86E-05|0.033171>,<7.31E-05|0.047349>,<9.77E-05|0.057941>,<0.000122|0.066835>,<0.000147|0.074711>,<0.000171|0.081864>,<0.000196|0.088463>,<0.000221|0.094621>,<0.000245|0.10042>,<0.00027|0.10591>,<0.000295|0.11114>,<0.000319|0.11614>,<0.000344|0.12091>,<0.000368|0.12544>,<0.000393|0.12977>,<0.000418|0.13391>,<0.000442|0.13788>,<0.000467|0.1417>,<0.000492|0.14539>,<0.000516|0.14896>,<0.000541|0.15244>,<0.000565|0.15582>,<0.00059|0.15911>,<0.000615|0.16231>,<0.000639|0.16544>,<0.000664|0.1685>,<0.000689|0.17148>,<0.000713|0.1744>,<0.000738|0.17725>,<0.000762|0.18005>,<0.000787|0.18279>,<0.000812|0.18548>,<0.000836|0.18811>,<0.000861|0.1907>,<0.000885|0.19324>,<0.00091|0.19574>,<0.000935|0.19819>,<0.000959|0.2006>,<0.000984|0.20297>,<0.001008|0.20531>,<0.001033|0.20761>,<0.001058|0.20987>,<0.001082|0.2121>,<0.001107|0.2143>,<0.001131|0.21647>,<0.001156|0.2186>,<0.00118|0.22071>,<0.001205|0.22279>,<0.00123|0.22484>,<0.001254|0.22687>,<0.001279|0.22886>,<0.001303|0.23083>,<0.001328|0.23278>,<0.001352|0.2347>,<0.001377|0.2366>,<0.001401|0.23847>,<0.001425|0.24032>,<0.00145|0.24215>,<0.001474|0.24396>,<0.001498|0.24575>,<0.001523|0.24751>,<0.001547|0.24926>,<0.001571|0.25099>,<0.001596|0.25277>,<0.00162|0.25433>,<0.001644|0.25571>,<0.001668|0.25708>,<0.001693|0.25843>,<0.001717|0.25978>,<0.001741|0.26112>,<0.001765|0.26245>,<0.001789|0.26377>,<0.001813|0.26508>,<0.001838|0.26638>,<0.001862|0.26767>,<0.001886|0.26895>,<0.00191|0.27022>,<0.001934|0.27148>,<0.001958|0.27273>,<0.001982|0.27396>,<0.002006|0.27519>,<0.00203|0.27641>,<0.002054|0.27762>,<0.002078|0.27881>,<0.002102|0.28>,<0.002126|0.28118>,<0.00215|0.28235>,<0.002174|0.28352>,<0.002198|0.28467>,<0.002222|0.28582>,<0.002246|0.28695>,<0.00227|0.28808>,<0.002294|0.2892>,<0.002318|0.29032>,<0.002342|0.29142>,<0.002366|0.29252>,<0.00239|0.29361>,<0.002414|0.29469>,<0.002438|0.29577>,<0.002462|0.29684>,<0.002486|0.2979>,<0.002509|0.29895>,<0.002533|0.3>,<0.002557|0.30104>,<0.002581|0.30207>,<0.002605|0.3031>,<0.002629|0.30412>,<0.002653|0.30513>,<0.002677|0.30614>,<0.002701|0.30714>,<0.002724|0.30813>,<0.002748|0.30912>,<0.002772|0.3101>,<0.002796|0.31107>,<0.00282|0.31204>,<0.002844|0.313>,<0.002868|0.31396>,<0.002891|0.31491>,<0.002915|0.31586>>:

RecenteredPowerFit:= proc(TW::Matrix, b_range::range(realcons))

local

     T:= TW[..,1], W:= TW[..,2],

     b:= Optimization:-NLPSolve(

          b-> Statistics:-PowerFit(T-~b, W, output= residualsumofsquares),

          b_range, method= branchandbound, objectivetarget= 0

     )[2][1],

     Res:= Statistics:-PowerFit(T-~b, W)

;

     ['A' = Res[1], ':-b' = b, 'c'= Res[2]]

end proc:

 

First try with the full data set. Smallest t value is 2e-5.

Sol_Carl1:= RecenteredPowerFit(TW, -1..TW[1,1] - 1e-7);

[A = HFloat(5.028381968959027), b = HFloat(1.9211938376704556e-5), c = HFloat(0.46614083419427993)]

Sol_Haley:= [A = 2.48, b = 1.1e-4, c = .35]:

PowerFitCheck:= proc(TW::Matrix, b::realcons)
local
     T:= TW[..,1] -~ b, W:= TW[..,2], n:= op([1,1],TW),
     nRSS:= Statistics:-PowerFit(T, W, output= residualsumofsquares)/n
;
     nRSS, Statistics:-PowerFit(T, W)
end proc:
     

PowerFitCheck(TW, eval(b, Sol_Carl1));

0.844786726509007e-3, Vector(2, {(1) = 5.02838196895903, (2) = .466140834194280})

Res:= PowerFitCheck(TW[5.., ..], eval(b, Sol_Haley));

Res := 0.215501674137240e-2, Vector(2, {(1) = 2.49659481877520, (2) = .354232604162157})

Sol_Haley5:= [A= Res[2][1], Sol_Haley[2], c=Res[2][2]];

[A = HFloat(2.496594818775202), b = 0.11e-3, c = HFloat(0.3542326041621568)]

Sol_Carl5:= RecenteredPowerFit(TW[5.., ..], -1..TW[5,1] - 1e-7);

[A = HFloat(3.682236898705571), b = HFloat(6.28308061404132e-5), c = HFloat(0.4158452235037357)]

PowerFitCheck(TW[5.., ..], eval(b, Sol_Carl5));

0.139855979031156e-3, Vector(2, {(1) = 3.68217591144777, (2) = .415842527560050})

P1:= plot(TW, style= point, symbol= cross, legend= data):

model:= A*(t-b)^c:

P:= plot([eval(model, Sol_Carl5), eval(model, Sol_Haley5)], t= 0..0.0029, legend= [Carl, Haley]):

plots:-display(P1,P);

 


Download RecenteredPowerFit.mw

@tomleslie The procedure that throws the error, `dsolve/numeric/DAE/make_proc`, is a monster: 792 lines and 36 parameters (as of Maple 18). The problem starts when floats in the equation are converted to exact rationals in line 22 and ends with the error message in line 719. I haven't traced what happens in between---it's a nightmare. It looks like the procedure has become a patchwork quilt over the years as features have been added to dsolve(..., numeric)

The same bug is manifested if the coefficient 3 is changed to 3.0.

First 483 484 485 486 487 488 489 Last Page 485 of 709