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

3 visualizaciones (últimos 30 días)
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 comentarios
Torsten
Torsten el 16 de Mayo de 2022
Editada: Torsten el 16 de Mayo de 2022
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
Roy Axford el 16 de Mayo de 2022
Editada: Matt J el 16 de Mayo de 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

Iniciar sesión para comentar.

Respuesta aceptada

Matt J
Matt J el 16 de Mayo de 2022
Editada: Matt J el 16 de Mayo de 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 comentario
Roy Axford
Roy Axford el 17 de Mayo de 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 :-)

Iniciar sesión para comentar.

Más respuestas (1)

Matt J
Matt J el 16 de Mayo de 2022
Editada: Matt J el 16 de Mayo de 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 comentario
Matt J
Matt J el 16 de Mayo de 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

Iniciar sesión para comentar.

Categorías

Más información sobre Matrices and Arrays en Help Center y File Exchange.

Productos


Versión

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by