Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
3  
try chisq[x_?NumericQ] !Mathematica graphics –  Nasser Jan 7 at 1:31
    
Thanks! Do you understand why it needs to be forced to use numerics in this case? –  Matt Reece Jan 7 at 4:05
    

2 Answers 2

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:

enter image description here

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]
share|improve this answer

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}]

enter image description here

{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.

share|improve this answer
    
Thanks. Using UnitStep[] is basically an improved version of the kluge I was doing with Tanh[]. (For some reason I was worried that a non-smooth function would trip it up.) I think Nasser's ?NumericQ trick is somewhat nicer, although I still wonder why Mathematica doesn't get this right from the outset. –  Matt Reece Jan 7 at 4:07
    
@MattReece It's because Mma tries to use chisq'[x] to find the minimum -- like every good calculus student ;) –  Michael E2 Jan 7 at 4:19

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.