Effacer les filtres
Effacer les filtres

caculate confidence interval from customized pdf

6 vues (au cours des 30 derniers jours)
苏越 徐
苏越 徐 le 19 Mar 2024
Commenté : 苏越 徐 le 19 Mar 2024
Hi
I'm wondering How can I caculate the confidence interval of customized pdf e.g. Gaussian mixture distribution?
pdf=@(x) w1*normpdf(x,mu1,sigma1)+w2*normpdf(x,mu2,sigma2);
cdf=@(x) integral(pdf,-Inf,x);
As icdf function only support specified distribution, I'm wondering how to caculate the shortest confidence interval?

Réponse acceptée

David Goodmanson
David Goodmanson le 19 Mar 2024
Modifié(e) : David Goodmanson le 19 Mar 2024
Hello SX,
Ordinarily to find the cdfs you would have to use numerical integration. In this case for the normal distributions, the cdf function is available. Then you can interpolate using the cdf as the independent variable. Here is an example. In the plot you get a wide minimum which you might expect.
mu1 = 1;
mu2 = 2;
sig1 = 1;
sig2 = 3;
w1 = .3;
w2 = .7;
c = .9; % confidence span, there is probably a better name for this
x = -20:.00001:20;
cdf = w1*normcdf(x,mu1,sig1) +w2*normcdf(x,mu2,sig2);
cdn = linspace(min(cdf),max(cdf)-c,1e4);
xdn = interp1(cdf,x,cdn);
cup = linspace(min(cdf)+c,max(cdf),1e4);
xup = interp1(cdf,x,cup);
figure(1); grid on
plot(xup-xdn)
[x0 ind] = min(xup-xdn);
xdn(ind) % lower end of confidence interval
xup(ind) % upper end of confidence interval
cdn(ind) % lower cdf value
cup(ind) % upper cdf value
% ans = -2.4087
% ans = 6.3858
% ans = 0.0497
% ans = 0.9497
D = xup(ind)-xdn(ind) % the result
cup(ind)-cdn(ind) % check, should be c = confidence span
% D = 8.7945
% ans = 0.9000
% try a different case, get a larger confidence interval
xtest = interp1(cdf,x,[.07 .97]);
Dtest = diff(xtest)
% xtest = -1.8602 7.1554
% Dtest = 9.0156
  3 commentaires
David Goodmanson
David Goodmanson le 19 Mar 2024
See what you think of the modified answer above
苏越 徐
苏越 徐 le 19 Mar 2024
Thank you David! It' good idea to narrow the range for ends of confidence span and then search for the shortest.

Connectez-vous pour commenter.

Plus de réponses (0)

Community Treasure Hunt

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

Start Hunting!

Translated by