ind2sub is not working properly.

I have the following piece of code where ind2sub is not working properly because when I plot the function PeLRTfusion in 3D, the results are not compatible. I want to find the coordinates of minimum value of PeLRTfusion. Also, when I use the resulted coordinates again in PeLRTfusion, the value of PeLRTfusion is different than the resulted min value from the code.

MuX1=1;
MuX2=2;
SigmaX1=1;
SigmaX2=0.5;
PH1=0.6;
PH0=1-PH1;
const=6;%A constant for how many sigmas away from the mean the program should go.
step=0.1;%To control the accuracy of the results.
Lower=-const*max(SigmaX1,SigmaX2);
Upper=const*max(SigmaX1,SigmaX2)+max(MuX1,MuX2);
LamdaSearchRange=Lower:step:Upper;
%Pe using the LRT fusion rule.
fun=@(Lamda1,Lamda2) PeLRTfusion(PH1,SigmaX1,SigmaX2,MuX1,MuX2,Lamda1,Lamda2);
PeOfLRT=bsxfun(fun,LamdaSearchRange,LamdaSearchRange');
[MinPeOfLRTval, MinIdx] = min(PeOfLRT(:));
MinPeOfLRTval
[MinLamda1Index, MinLamda2Index] = ind2sub(size(PeOfLRT),MinIdx);
MinLamda1=LamdaSearchRange(MinLamda1Index)
MinLamda2=LamdaSearchRange(MinLamda2Index)

Any help is really appreciated. Thank you.

[Image Analyst]

Since I can't test it because you forgot to include PeLRTfusion, I gave up trying to help. Include PeLRTfusion if you want people to be able to run your code and fix it.

[Hadi Kasasbeh]

Sorry for that, here is the code

function Pe=PeLRTfusion(PH1,SigmaX1,SigmaX2,MuX1,MuX2,Lamda1,Lamda2)
PH0=1-PH1;
lamda=PH0/PH1;
PX1e1H1=qfunc((Lamda1-MuX1)/SigmaX1);%P(X_1^q=1|H1)
PX1e0H1=1-PX1e1H1;%P(X_1^q=0|H1)
PX2e1H1=qfunc((Lamda2-MuX2)/SigmaX2);%P(X_2^q=1|H1)
PX2e0H1=1-PX2e1H1;%P(X_2^q=0|H1)
PX1e1H0=qfunc(Lamda1/SigmaX1);%P(X_1^q=1|H0)
PX1e0H0=1-PX1e1H0;%P(X_1^q=0|H0)
PX2e1H0=qfunc(Lamda2/SigmaX2);%P(X_2^q=1|H0)
PX2e0H0=1-PX2e1H0;%P(X_2^q=0|H0)
PD1=bi2de((PX1e0H1*PX2e0H1./(PX1e0H0*PX2e0H0))>=lamda);%P(D=1|X_1^q=0,X_2^q=0)
PD2=bi2de((PX1e0H1*PX2e1H1./(PX1e0H0*PX2e1H0))>=lamda);%P(D=1|X_1^q=0,X_2^q=1)
PD3=bi2de((PX1e1H1*PX2e0H1./(PX1e1H0*PX2e0H0))>=lamda);%P(D=1|X_1^q=1,X_2^q=0)
PD4=bi2de((PX1e1H1*PX2e1H1./(PX1e1H0*PX2e1H0))>=lamda);%P(D=1|X_1^q=1,X_2^q=1)
Pe=PH1+PD1.*(PH0*PX1e0H0*PX2e0H0-PH1*PX1e0H1*PX2e0H1)+...
    PD2.*(PH0*PX1e0H0*PX2e1H0-PH1*PX1e0H1*PX2e1H1)+...
    PD3.*(PH0*PX1e1H0*PX2e0H0-PH1*PX1e1H1*PX2e0H1)+...
    PD4.*(PH0*PX1e1H0*PX2e1H0-PH1*PX1e1H1*PX2e1H1);

Thank you.

[Hadi Kasasbeh]

This the 3D plotting part of the first code

fsurf(@(Lamda1,Lamda2) PeLRTfusion(PH1,SigmaX1,SigmaX2,MuX1,MuX2,Lamda1,Lamda2),[(-const*SigmaX1) (max(MuX1,MuX2)+const*SigmaX1)])
xlabel('\lambda_1'); ylabel('\lambda_2'); zlabel('P_e');

[Stephen23]

ind2sub seems to be working perfectly:

[MinPeOfLRTval, MinIdx] = min(PeOfLRT(:))
MinPeOfLRTval =
 -1.6092e+042
MinIdx =
        4018
>> val = PeOfLRT(MinLamda1Index,MinLamda2Index)
val =
 -1.6092e+042
>> MinPeOfLRTval - val
ans =
     0

[Hadi Kasasbeh]

Thank you Stephen for your reply, but this is the output of my program:

MinPeOfLRTval = 0.0227

MinLamda1Index =70

MinLamda2Index = 29

MinLamda1 = 0.9000

MinLamda2 = -3.2000

Now, if you use Lamda1=MinLamda1 = 0.9000 and Lamda2=MinLamda2 = -3.2000 in the function PeOfLRT, you get the value 0.3497 which in not the same as 0.0227. Also if you look to the plot of PeOfLRT you can see clearly that the minimum point is Not (0.9,-3.2).

Answer 1

[Stephen23]

Your function PeLRTfusion returns different outputs values depending on the size of the inputs:

>> [v,x] = min(fun(LamdaSearchRange(29),LamdaSearchRange)) % scalar & vector
v =
 -1.6092e+042
x =
    70
>> [v,x] = min(fun(LamdaSearchRange(29),LamdaSearchRange(70))) % scalar & scalar
v =
     0.022712
x =
     1

So when you try it with two scalars it gives a different output than with a scalar and a vector. My installed help for bsxfun states "if an M-file function is specified, it must be able to accept either two column vectors of the same size, or one column vector and one scalar, and return as output a column vector of the size as the input values", which tell us that MATLAB passes the entire columns at once to fun. Unfortunately they seem to have removed this useful info from the online help.

In any case, you could either:

1. fix your function to work correctly when called with column lambdas.
2. use arrayfun (which works element-by-element), instead of bsxfun (which does not).

For the second option you simply need to replace this:

PeOfLRT=bsxfun(fun,LamdaSearchRange,LamdaSearchRange');

with this:

[L1,L2] = ndgrid(LamdaSearchRange);
PeOfLRT = arrayfun(fun,L1,L2);

which gives these values:

MinPeOfLRTval =
     0.022712
MinIdx =
        9758
MinLamda1Index =
    29
MinLamda2Index =
    70

We can check these values by hand:

>> fun(LamdaSearchRange(29),LamdaSearchRange(70))
ans =
     0.022712

and by viewing the values in 3D:

>> surf(PeOfLRT)

Judging by the long trough there many be many values that are very close to the minimum. Lets check:

>> nnz(abs(PeOfLRT-min(PeOfLRT(:)))<1e-5)
ans =
   141
>> nnz(abs(PeOfLRT-min(PeOfLRT(:)))<1e-6)
ans =
     6

Which means that quite possibly the minimum value found depends on some floating point error more than being the "exact" function minimum.

[Hadi Kasasbeh]

It worked perfectly! Thank you Stephen so much. I also guess the title of this question should be changed accordingly.