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

3 views (last 30 days)
Roy Axford
Roy Axford on 16 May 2022
Commented: Roy Axford on 17 May 2022
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?
  2 Comments
Roy Axford
Roy Axford on 16 May 2022
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

Sign in to comment.

Accepted Answer

Matt J
Matt J on 16 May 2022
Edited: Matt J on 16 May 2022
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
  1 Comment
Roy Axford
Roy Axford on 17 May 2022
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 :-)

Sign in to comment.

More Answers (1)

Matt J
Matt J on 16 May 2022
Edited: Matt J on 16 May 2022
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)
  1 Comment
Matt J
Matt J on 16 May 2022
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

Sign in to comment.

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!