What is the best way set the search interval used by fminbnd?

I am using fminbnd to find the minimum of a function in which the independent variable is a probability. Therefore, I thought I could use 0 and 1 as the bounds of the search interval. However, this approach is unrelaible and I have had to "cheat" by setting the search interval to something that more tightly bounds the known correct answers for my test cases. Thus, I have not built confidence in my ability to use fminbnd to find correct answers that I don't already know.

Have I missed the guidance for setting the search interval?

[Torsten]

We must see your application to comment on this.

If you have a discrete probability distribution, only discrete values for probability are taken between 0 and 1. Thus a continuous minimizer like fminbnd is not applicable.

[Roy Axford]

Torsten,

I am dealing with the binomial distribution. The probabilites are not constrained to discrete values whereas the numbers of trials and successes are. So, I'm not sure that using fminbnd is invalid. It does give correct answers if I sufficiently constrain the search region, but I want to get away from having to do that as I said in my original post.

Here are my two files, a function that uses binopdf and a script that uses fminbnd.

Thanks for any advice you can provide.

Regards,

Roy

=== BEGIN SCRIPT ===

% test_summing_binomial_pdf_2
% Roy Axford
% May 16, 2022
%
% Want to use fminunc to find the value of p that minimizes (finds the zero of) sum_prob_minus_epsilon_half.
%
% Set number of trials and errors, and confidence coefficient.
n = 1000;     % Number of trials.  Use 1000.
S = 1;        % Number of errors.  Use 1, 10, & 100.
beta = 0.95;  % Confidence coefficient.
epsilon = 1-beta;
epsilon_half = epsilon/2;
% p1 SECTION - Lower end of confidence interval.
firstindex = S;
lastindex = n;
x_vector = firstindex:lastindex;
x_vector = x_vector';
f = @(p)sumbinopdf_equals_epsilon_half_2(p,x_vector,n,epsilon_half);
% Use optimset to set options used by the solver 'fminbnd'.
options = optimset('MaxFunEvals',10000,'MaxIter',10000,'TolX',1.0e-8);
p1 = fminbnd(f,1e-5,5e-5,options)
% 
% p2 SECTION - Upper end of confidence interval.
firstindex = 0;
lastindex = S;
x_vector = firstindex:lastindex;
x_vector = x_vector';
f = @(p)sumbinopdf_equals_epsilon_half_2(p,x_vector,n,epsilon_half);
p2 = fminbnd(f,3e-3,7e-3,options)
disp('Script finished.')

=== END SCRIPT ===

function [the_difference] = sumbinopdf_equals_epsilon_half_2(p,x_vec,n,epsilon_half)
% Sums the binomial pdf, binopdf(x,n,p) over x_vec and compare the sum to epsilon/2.
% Intent: Use this function to search for a value of p that satisfies
% sum(binopdf(x_vec,n_vec,p_vec)) = epsilon/2
% in the context of finding confidence limits p1 & p2 as decribed in
% Section 14.54 of Burington & May, Handbook of Probability and Statistics
% with Tables, 1958.
% Formulate the difference between sum(binopdf(x_vec,n_vec,p_vec)) &
% epsilon_half so the desired result is that this function equals zero.
%
% Inputs
% x_vec:        vector, a range of number of successes.
% n:            scalar, number of trials.
% p:            scalar, p=Pr[success] - script varies this in calls to this function.
% epsilon_half: scalar = 0.5*(1-beta); beta is the confidence coefficient
%
% Output
% the_difference: absolute value of the difference between the desired sum of probabilities and epsilon_half.
%
n_vec = n.*ones(size(x_vec));
p_vec = p.*ones(size(x_vec));
the_difference = abs(sum(binopdf(x_vec,n_vec,p_vec))-epsilon_half);
end

Answer 1

[Matt J]

Your minimization problem is really a root-finding problem in disguise. It is better to use fzero for such things. As you can see below, we get good results over a wider search interval:

n = 1000;     % Number of trials.  Use 1000.
S = 1;        % Number of errors.  Use 1, 10, & 100.
beta = 0.95;  % Confidence coefficient.
epsilon = 1-beta;
epsilon_half = epsilon/2;
% p1 SECTION - Lower end of confidence interval.
firstindex = S;
lastindex = n;
x_vector = firstindex:lastindex;
x_vector = x_vector';
f = @(p)sumbinopdf_equals_epsilon_half_2(p,x_vector,n,epsilon_half);

[p1,fval,exitflag]=fzero(f,[0,1])

p1 = 2.5317e-05

fval = -6.1944e-14

exitflag = 1

% p2 SECTION - Upper end of confidence interval.
firstindex = 0;
lastindex = S;
x_vector = firstindex:lastindex;
x_vector = x_vector';
f = @(p)sumbinopdf_equals_epsilon_half_2(p,x_vector,n,epsilon_half);

[p2,fval,exitflag]=fzero(f,[0,1])

p2 = 0.0056

fval = 1.0408e-16

exitflag = 1

function [the_difference] = sumbinopdf_equals_epsilon_half_2(p,x_vec,n,epsilon_half)
% Sums the binomial pdf, binopdf(x,n,p) over x_vec and compare the sum to epsilon/2.
% Intent: Use this function to search for a value of p that satisfies
% sum(binopdf(x_vec,n_vec,p_vec)) = epsilon/2
% in the context of finding confidence limits p1 & p2 as decribed in
% Section 14.54 of Burington & May, Handbook of Probability and Statistics
% with Tables, 1958.
% Formulate the difference between sum(binopdf(x_vec,n_vec,p_vec)) &
% epsilon_half so the desired result is that this function equals zero.
%
% Inputs
% x_vec:        vector, a range of number of successes.
% n:            scalar, number of trials.
% p:            scalar, p=Pr[success] - script varies this in calls to this function.
% epsilon_half: scalar = 0.5*(1-beta); beta is the confidence coefficient
%
% Output
% the_difference: absolute value of the difference between the desired sum of probabilities and epsilon_half.
%
n_vec = n.*ones(size(x_vec));
p_vec = p.*ones(size(x_vec));
the_difference = sum(binopdf(x_vec,n_vec,p_vec))-epsilon_half;
end

[Roy Axford]

Matt J.,

Thanks very much. Yes, fzero works much better for this application. I hadn't used fzero before and I was unaware of it.

I now have very solid code for solving this problem, and I expect to find more uses for fzero in the future.

All the best,

Roy :-)

Answer 2

[Matt J]

and I have had to "cheat" by setting the search interval to something that more tightly bounds the known correct answers for my test cases.

In cases where you don't know the solution, you can sweep the interval to find a tighter subinterval, like in the following:

xsamples=linspace(0,1,10);
fsamples = arrayfun(fun, xsamples);

[~,im]=min(fsamples);

a=xsamples(  max( 1,im-1)   );
b=xsamples(  min(10,im+1)   );

x=fminbnd(fun, a,b)

[Matt J]

Applying this to your problem seems to do alright:

n = 1000;     % Number of trials.  Use 1000.
S = 1;        % Number of errors.  Use 1, 10, & 100.
beta = 0.95;  % Confidence coefficient.
epsilon = 1-beta;
epsilon_half = epsilon/2;
% p1 SECTION - Lower end of confidence interval.
firstindex = S;
lastindex = n;
x_vector = firstindex:lastindex;
x_vector = x_vector';
f = @(p)sumbinopdf_equals_epsilon_half_2(p,x_vector,n,epsilon_half);

[p1,fval,exitflag] = searchRoot(f,30)

p1 = 2.5317e-05

fval = 2.7794e-10

exitflag = 1

firstindex = 0;
lastindex = S;
x_vector = firstindex:lastindex;
x_vector = x_vector';
f = @(p)sumbinopdf_equals_epsilon_half_2(p,x_vector,n,epsilon_half);


[p2,fval,exitflag] = searchRoot(f,30)

p2 = 0.0056

fval = 3.5502e-10

exitflag = 1

function varargout = searchRoot(fun,N)

xsamples=linspace(0,1,N);
fsamples = arrayfun(fun, xsamples);

[~,im]=min(fsamples);

a=xsamples(  max( 1,im-1)   );
b=xsamples(  min(N,im+1)   );

[varargout{1:nargout}]=fminbnd(fun, a,b,...
                       optimset('TolX',1e-12,'MaxIter',1e20,'MaxFunEvals',1e20));

end

function [the_difference] = sumbinopdf_equals_epsilon_half_2(p,x_vec,n,epsilon_half)
n_vec = n.*ones(size(x_vec));
p_vec = p.*ones(size(x_vec));
the_difference = abs(sum(binopdf(x_vec,n_vec,p_vec))-epsilon_half);
end