FindMinimum gets confused by If[...]

Question

I'm confused by why the following doesn't work. I'm trying to numerically minimize a function that has an If statement in its definition, as follows:

means = {3.0, 1.5, 2.1};    
sigmashi = {0.5, 0.4, 0.3};
sigmaslo = {0.4, 0.6, 0.25};
coeffA = {1.2, 3.2, 2.0};
coeffB = {0.3, -0.4, 0.5};
chisq[x_] := 
 Sum[If[coeffA[[i]] x + coeffB[[i]] > 
    means[[i]], (coeffA[[i]] x + coeffB[[i]] - means[[i]])^2/
    sigmashi[[i]]^2, (coeffA[[i]] x + coeffB[[i]] - means[[i]])^2/
    sigmaslo[[i]]^2], {i, 1, Length[means]}]    
FindMinimum[chisq[x], {x, 1.0}]

This returns a long string of error messages:

Part::pspec: Part specification i is neither a machine-sized integer nor a list of machine-sized integers. >> Part::pspec: Part specification i is neither a machine-sized integer nor a list of machine-sized integers. >> Part::pspec: Part specification i is neither a machine-sized integer nor a list of machine-sized integers. >> General::stop: Further output of Part::pspec will be suppressed during this calculation. >> FindMinimum::nrnum: The function value ({0.3,-0.4,0.5}[[i]]+1. <<1>>-{3.,1.5,2.1}[[i]])^2/{0.4,0.6,0.25}[[i]]^2+(2 ({0.3,-0.4,0.5}[[i]]+<<1>>-{3.,1.5,2.1}[[i]])^2)/{0.5,0.4,0.3}[[i]]^2 is not a real number at {x} = {1.}. >>

Why doesn't FindMinimum like my function? Why does it even care about the "i"? Is there a way to force Mathematica to do the minimization purely numerically, and not care about how the function computes its result? (NMinimize returns the same errors.)

I can kluge my way around this with the following alternative:

myif[z_, x_, y_] := ((Tanh[100 z] + 1)/2) x + ((Tanh[-100 z] + 1)/2) y
chisqA[x_] := 
 Sum[myif[coeffA[[i]] x + coeffB[[i]] - 
    means[[i]], (coeffA[[i]] x + coeffB[[i]] - means[[i]])^2/
    sigmashi[[i]]^2, (coeffA[[i]] x + coeffB[[i]] - means[[i]])^2/
    sigmaslo[[i]]^2], {i, 1, Length[means]}]

But it isn't pretty.

Answer 1

Using the initiative of bills and rewriting your function:

csq[x_] := 
 Total[UnitStep[x #1 + #2 - #3] (#1 x + #2 - #3)^2/#4^2 + 
     UnitStep[-(x #1 + #2) + #3] (#1 x + #2 - #3)^2/#5^2 & @@@ 
   Transpose[{coeffA, coeffB, means, sigmashi, sigmaslo}]]

then:

NMinimize[csq[x], x]

yields:

{21.6448, {x -> 0.798905}}

Or

FindMinimum[csq[x], {x, 1}]

yields:

{21.6448, {x -> 0.798905}}

Just for confirmation: plotting original function and min, redefined function and min, difference between two functions:

Code:

min = FindMinimum[csq[x], {x, 1}];
GraphicsRow[{Plot[chisq[x], {x, -3, 3}, PlotLabel -> "chisq(x)", 
   Epilog -> {Red, PointSize[0.04], 
     Point[{x /. min[[2]], min[[1]]}]}], 
  Plot[csq[x], {x, -3, 3}, 
   Epilog -> {Red, PointSize[0.04], Point[{x /. min[[2]], min[[1]]}]},
    PlotLabel -> "csq(x)"], 
  Plot[chisq[x] - csq[x], {x, -3, 3}, 
   PlotLabel -> "chisq(x)-csq[x]"]}, ImageSize -> 700]

Answer 2

Let's see if we can figure out what's really causing the problem. Here's a simplified version of your problem:

f[x_] := If[x > 0, x^2, x^4];
Plot[f[x], {x, -1, 1}]
FindMinimum[f[x], {x, 0.1}]

{0., {x -> 0.}}

So there's no problem really If alone. Let's try with the indices into the lists:

a = {1, 2}; b = {3, 4};
g[x_] := Sum[If[a[[i]] x > 0, a[[i]] x^2, b[[i]] x^4], {i, 1, 2}];
FindMinimum[g[x], {x, 0.1}]
FindMinimum::nrnum: "The function value 0.02\ {1,2}[[i]] is not a real number at {x} = {0.1`}"

So there is the problem. Here's a way to fix it, using the UnitStep function:

a = {1, 2}; b = {3, 4};
h[x_] := Sum[UnitStep[a[[i]] x]*a[[i]] x^2 + UnitStep[-a[[i]] x]*b[[i]] x^4, {i, 1, 2}];
FindMinimum[h[x], {x, 0.1}]
{2.23889*10^-29, {x -> -4.22896*10^-8}}

I'm sure you can figure out how to go back and apply UnitStep to your function.