Actually, I'd suggest using Fit. The expression you're fitting is linear in the parameters, so LinearFit can handle it; and LinearFit finds the globally optimal parameter values, as opposed to NonlinearFit which finds only a locally optimal point. The Fit command selects LinearFit if appropriate and NonlinearFit if necessary, and I'd recommend using that when you're not absolutely certain that you need NonlinearFit. That would involve just replacing the command NonlinearFit with Fit.

In this case it actually doesn't make a difference for the values for a and b that you get out, by the way, but it's just a more reliable process this way.

Oh, and if you need to extract the columns of an Array:

V := Data[.., 1];
C := Data[.., 2];

Hope this helps,

Erik Postma
Maplesoft.

Actually, I'd suggest using Fit. The expression you're fitting is linear in the parameters, so LinearFit can handle it; and LinearFit finds the globally optimal parameter values, as opposed to NonlinearFit which finds only a locally optimal point. The Fit command selects LinearFit if appropriate and NonlinearFit if necessary, and I'd recommend using that when you're not absolutely certain that you need NonlinearFit. That would involve just replacing the command NonlinearFit with Fit.

In this case it actually doesn't make a difference for the values for a and b that you get out, by the way, but it's just a more reliable process this way.

Oh, and if you need to extract the columns of an Array:

V := Data[.., 1];
C := Data[.., 2];

Hope this helps,

Erik Postma
Maplesoft.

@pagan : This is amazing. Somehow, Firefox (assisted by some scripting trickery, for sure) must have thought it would be a great idea to use that pasted expression, write it to my temporary directory on the local filesystem as /tmp/moz-screenshot.png, and then insert the img tag referring to the local uri file:///tmp/moz-screenshot.png. Wow. I just checked and indeed I have that file on my file system here.

Of course I had no idea that that's what's going on. I just saw an image appear in my post... So anyway, it looks like this is not very useful.

Erik.

@pagan : This is amazing. Somehow, Firefox (assisted by some scripting trickery, for sure) must have thought it would be a great idea to use that pasted expression, write it to my temporary directory on the local filesystem as /tmp/moz-screenshot.png, and then insert the img tag referring to the local uri file:///tmp/moz-screenshot.png. Wow. I just checked and indeed I have that file on my file system here.

Of course I had no idea that that's what's going on. I just saw an image appear in my post... So anyway, it looks like this is not very useful.

Erik.

I tried a few random (positive) values for the other parameters, and for a couple of them, there was a unique real solution for i, but not with these values:

y = 2
a = 9
b = 1/10
c = 1/10
d = 5
g = 3/4
h = 1/8

Here, there are still a number of solutions (Maple finds 51 of them), but they're all complex, and I guess that would not be valid in the context of neural networks. Here's a plot of the right hand side of the equation minus y:

And indeed, we can see that if y > g + h, then substituting anything positive for the exponentials is not going to lead to a solution. So this seems to indicate that the function is not invertible -- unless this violates some constraint on the values that the parameters should have?

Erik.

I tried a few random (positive) values for the other parameters, and for a couple of them, there was a unique real solution for i, but not with these values:

y = 2
a = 9
b = 1/10
c = 1/10
d = 5
g = 3/4
h = 1/8

Here, there are still a number of solutions (Maple finds 51 of them), but they're all complex, and I guess that would not be valid in the context of neural networks. Here's a plot of the right hand side of the equation minus y:

And indeed, we can see that if y > g + h, then substituting anything positive for the exponentials is not going to lead to a solution. So this seems to indicate that the function is not invertible -- unless this violates some constraint on the values that the parameters should have?

Erik.

@Alec Mihailovs: yeah, that's my doing. I had to fiddle around a bit to get the command to show up right (for certain values of right - the 0 and 3 are in an italic font, not roman, but I'd consider that minor). I guess there are some kinks to be worked out of the code dealing with that. For the record, here's the intended command in monospaced font:

plot(sin(t^2), t = 0 .. 3, labels = [sqrt(time), amplitude]);

Thanks, Erik.

@Alec Mihailovs: yeah, that's my doing. I had to fiddle around a bit to get the command to show up right (for certain values of right - the 0 and 3 are in an italic font, not roman, but I'd consider that minor). I guess there are some kinks to be worked out of the code dealing with that. For the record, here's the intended command in monospaced font:

plot(sin(t^2), t = 0 .. 3, labels = [sqrt(time), amplitude]);

Thanks, Erik.

@Alec Mihailovs : that's right, and that's a good point to make. I just wanted to make the slightly different point that it might not be a great idea to do RREF at all.

My 1/3 example was bad; what I meant is that such fractions can come up during the computation, not just if they're present in the input (or even in the correct output): if you encounter a pivot during RREF that you cannot divide exactly by in the arithmetic that's being used (e.g. 3 as a pivot for either binary or decimal arithmetic), then other entries in that row which are not a mulitple of that pivot will be divided by that pivot and will not be exact anymore. An example would be this matrix:

Actually Maple finds the right RREF here, but if you do it naively, then scaling the first row leaves you with <1., 1.3333333332, -1.6666666665> (in the binary arithmetic) and you get the wrong answer in the end, even though both input and (correct) output contain only things that are correctly representable.

Erik Postma
Maplesoft.

@Alec Mihailovs : that's right, and that's a good point to make. I just wanted to make the slightly different point that it might not be a great idea to do RREF at all.

My 1/3 example was bad; what I meant is that such fractions can come up during the computation, not just if they're present in the input (or even in the correct output): if you encounter a pivot during RREF that you cannot divide exactly by in the arithmetic that's being used (e.g. 3 as a pivot for either binary or decimal arithmetic), then other entries in that row which are not a mulitple of that pivot will be divided by that pivot and will not be exact anymore. An example would be this matrix:

Actually Maple finds the right RREF here, but if you do it naively, then scaling the first row leaves you with <1., 1.3333333332, -1.6666666665> (in the binary arithmetic) and you get the wrong answer in the end, even though both input and (correct) output contain only things that are correctly representable.

Erik Postma
Maplesoft.

It turns out that most of the time is spent in overhead in the Statistics package. There is fast external C code for shuffling, but finding that code takes a while. If you use the debugger to step through a Statistics[Shuffle] call, you can see that in the end, it calls an external function called MapleShuffle. You can shave off another factor of 3 or so by calling that directly. The following example was directly adapted from your code; it may of course not work in past or future versions of Maple (it does work in Maple 14, at least).

restart;
with(Statistics):
opacity := kernelopts(opaquemodules=false):
Statistics:-ExternalSupport:-initializeHwLibrary():
shuffle := Statistics:-ExternalSupport:-DefineExternal('MapleShuffle');
kernelopts(opaquemodules=opacity):
interface(rtablesize=12):
L := Array([2, 2, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0]):
Y := Array(1..6,[seq(0,i=2..7)]):
X := Array(1..6,[seq(0,i=2..7)]):
n:=792000:
# Digits:=2: # why?
st:=time():
for i to n do
L := shuffle(L):
E1:=add(L[j],j=1..5):
X[E1-1] := X[E1-1]+1:
E2:=add(L[j],j=6..10):
Y[E1-1] := Y[E1-1] + E2:
end do:
time()-st;

Hope this helps,

Erik Postma
Maplesoft.

It turns out that most of the time is spent in overhead in the Statistics package. There is fast external C code for shuffling, but finding that code takes a while. If you use the debugger to step through a Statistics[Shuffle] call, you can see that in the end, it calls an external function called MapleShuffle. You can shave off another factor of 3 or so by calling that directly. The following example was directly adapted from your code; it may of course not work in past or future versions of Maple (it does work in Maple 14, at least).

restart;
with(Statistics):
opacity := kernelopts(opaquemodules=false):
Statistics:-ExternalSupport:-initializeHwLibrary():
shuffle := Statistics:-ExternalSupport:-DefineExternal('MapleShuffle');
kernelopts(opaquemodules=opacity):
interface(rtablesize=12):
L := Array([2, 2, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0]):
Y := Array(1..6,[seq(0,i=2..7)]):
X := Array(1..6,[seq(0,i=2..7)]):
n:=792000:
# Digits:=2: # why?
st:=time():
for i to n do
L := shuffle(L):
E1:=add(L[j],j=1..5):
X[E1-1] := X[E1-1]+1:
E2:=add(L[j],j=6..10):
Y[E1-1] := Y[E1-1] + E2:
end do:
time()-st;

Hope this helps,

Erik Postma
Maplesoft.

Thanks Alec. I always do that, and it always shows up perfectly for me in the preview, using either the HTML options. Then in the posted comment itself, boom, it's gone. Maybe there's a gremlin in the system that picks on my browser/OS combination or something.

Erik.

Thanks Alec. I always do that, and it always shows up perfectly for me in the preview, using either the HTML options. Then in the posted comment itself, boom, it's gone. Maybe there's a gremlin in the system that picks on my browser/OS combination or something.

Erik.

Here's one option, that works especially pleasantly with Acer's evalhf'able integrand (thanks Acer) :

n := 7;
plot(numericIntegral, 3 .. 7);

Now it may look like there's a zero at t = 5, but if you make Maple look closely, you can see that there's a reason why fsolve won't tell you that:

plot(numericIntegral, 0 .. 5);

It shows you that the integral is positive but extremely small, probably all due to roundoff error.

OK - I can't help myself, I have to give you an option to do it non-graphically :) If you're interested in just finding a spot where the integral is close to zero, you might want to use some custom code (instead of fsolve or RootFinding[NextZero]). For example, here's a very simple piece of code (using bisection) that finds the maximum value (up to an x tolerance, xtol) in a range where the supplied function is within a tolerance (ytol) of zero - or at least, a value that is a right end point of a range where the supplied function is near zero, there could be other such ranges further to the right:

maxApproximateRoot := proc(f :: procedure, r :: range(numeric), xtol :: numeric, ytol :: numeric, $)
local a, b, c, fa, fb, fc;
  a, b := op(evalf(r));
  fb := f(b);
  if abs(fb) <= ytol then
    return b;
  end if;

  fa := f(a);
  if fa * fb > 0 and abs(fa) > ytol then
    error "f has the same sign on both sides of the supplied interval and is not close to zero on either end - aborting";
  end if;

  # Invariant: (fa * fb < 0 or abs(fa) <= ytol) and abs(fb) > ytol
  while b - a > xtol do
    c := (a + b)/2;
    fc := f(c);

    if abs(fc) <= ytol or fc * fb < 0 then
      a := c;
      fa := fc;
    else
      b := c;
      fb := fc;
    end if;
  end do;

  return (a + b)/2;
end proc;

You could use it as follows:

maxApproximateRoot(numericIntegral, 3 .. 7, 1e-5, 1e-5);
# returns 5.019527435 (for n = 7)

Indeed, running

plot(numericIntegral, 5.01 .. 5.03);

gives a very interesting, non-continuous plot... numeric issues again!

Hope this helps,

Erik Postma
Maplesoft.