# Invalid computation of maxima

I am finding the following minimum

FindMinimum[{
Abs[
ArgMax[f[q, n, b, b0, b1, t, s], b] - b1
]
}, {b1, -0.5}]


given q=3/5, b0=1, t=3/2, s=1 and obtain the following graph for n = {5,...,25}:

My question is about the values at 13, 16, and 24. Why am I having these irregular and apparently wrong values and is there a way to get rid of this problem? Is this a problem of precision? Also, is it possible to make this computation faster?

Update Upon requests in comments, I give the definitions here. Function f is

f[q_, n_, b_, b0_, b1_, t_, s_] := Sum[(Binomial[n - 1, k]*t^(n - 1 - k)/(1 + t)^(n - 1)*g[q, n, b, k*b0 + (n - 1 - k)*b1, s]), {k, 0, n - 1}]


where function g is

g[q_, n_, b_, B_, s_] = (q*(Probability[x + b + B + 1/2 s^2 Log[q/(1 - q)] > 0,
x \[Distributed] NormalDistribution[n, s*n^(1/2)]]) + (1 -
q)*(Probability[x + b + B + 1/2 s^2 Log[q/(1 - q)] < 0,
x \[Distributed] NormalDistribution[-n, s*n^(1/2)]]))


Here is the plotting procedure, to make it easy for your replication. I set the function:

bb[q_, n_, b0_, t_, s_]:=b1/.Last[FindMinimum[{Abs[ArgMax[f[q,n,b,b0,b1,t,s],b]-b1]},{b1,-0.5}]]


and then plot with

DiscretePlot[bb[3/5,m,1,3/2,1],{m,5,25}]


In place of FindMinimum I used FindRoot as in

FindRoot[ArgMax[f[q,n,b,b0,b1,t,s],b]-b1==0,{b1,-0.5}]


and also, in another try, I got replaced Abs with conditions as in

FindMinimum[{ArgMax[f[q,n,b,b0,b1,t,s],b]-b1,ArgMax[f[q,n,b,b0,b1,t,s],b]-b1>0},{b1,-0.5}]


both of which did not help.

Finally, to show another proof that the problem with this deterministic computation is totally random, I present here one more plot, where I replace 1 with 1.0001 for parameter s.

Here is when s=1.0001

• Could you explain "an expected normal probability of binomial values" further? – Dr. belisarius Feb 16 '16 at 15:06
• @Dr.belisarius It is actually Sum[(Binomial[n - 1, k]*t^(n - 1 - k)/(1 + t)^(n - 1)*g[q, n, b, k*b0 + (n - 1 - k)*b1, s]), {k, 0, n - 1}] where the function g is really an error function (normal probability). – juror Feb 16 '16 at 15:13
• edit the question to include the definition of f and g – george2079 Feb 16 '16 at 16:22
• Please don't change the tags ... again – Dr. belisarius Feb 16 '16 at 16:27
• @Dr.belisarius Ok. – juror Feb 16 '16 at 16:29

There are a few ways to improve on this situation. SOme involve making these functions into "black boxes" that only work for explicit numeric input. That removes all symbolic preprocessing, along with warning messages that get produced there. Another is to force use of a the interior point method. This can be done by adding a constraint e.g. b1<=0. Or use Newton's method; that also seems to give a viable result.

For either of the above there is an important subtlety. They like differentiable functions. We have an Abs defining the objective function. Since it is a discrepancy we are measuring, we can just as well use the square of the value. In here lurks another subtlety: we want the squaring to be explicitly seen in FindMinimum and not part of the black-box function.

Also it seems like a good idea to use the numeric arg-max since it will be faster. The reson is that "exact" ArgMax will punt to NMinimize when given approximate input, whereas NArgMax will invoke the typically faster FindMaximum.

So here is the code.

q = 3/5;
b0 = 1;
t = 3/2;
s = 1.0001;

f[q_, n_?NumberQ, b_, b0_, b1_, t_, s_] :=
Sum[(Binomial[n - 1, k]*t^(n - 1 - k)/(1 + t)^(n - 1)*
g1[q, n, b, k*b0 + (n - 1 - k)*b1, s]), {k, 0, n - 1}]

g[q_, n_, b_, B_,
s_] = (q*(Probability[x + b + B + 1/2 s^2 Log[q/(1 - q)] > 0,
x \[Distributed] NormalDistribution[n, s*n^(1/2)]]) + (1 -
q)*(Probability[x + b + B + 1/2 s^2 Log[q/(1 - q)] < 0,
x \[Distributed] NormalDistribution[-n, s*n^(1/2)]]));

g1[q_, n_?NumberQ, b_, B_, s_] := g[q, n, b, B, s]

argmaxfN[q_, n_?NumberQ, b_, b0_, b1_?NumberQ, t_,
s_] := (NArgMax[f[q, n, b, b0, b1, t, s], b] - b1)


So here goes.

AbsoluteTiming[
res = Table[
FindMinimum[argmaxfN[q, n, b, b0, b1, t, s]^2, {b1, -0.5},
Method -> "Newton"], {n, 5, 25}];]

<a few warnings snipped>

(* Out[108]= {53.415487, Null} *)

ListPlot[b1 /. res[[All, 2]]]


As noted earlier, one can force usage of interior point code by adding a constraint.

AbsoluteTiming[
res = Table[
FindMinimum[{argmaxfN[q, n, b, b0, b1, t, s]^2,
b1 <= 0}, {b1, -0.5}], {n, 5, 25}];]

During evaluation of In[110]:= NArgMax::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations. >>

During evaluation of In[110]:= NArgMax::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations. >>

During evaluation of In[110]:= NArgMax::cvmit: Failed to converge to the requested accuracy or precision within 100 iterations. >>

During evaluation of In[110]:= General::stop: Further output of NArgMax::cvmit will be suppressed during this calculation. >>

(* Out[110]= {54.747887, Null} *)


The plot is pretty much the same as before, so I omit it.