Getting 1 sigma from non normal distribution

Hi all,

I'm trying to find the one sigma limits for a particular probability distribution (see image). I don't have a functional form of the distribution-it's generated from discrete data points-although I think I could find a reasonable match using cftool.

What I would like to do is find the 1 sigma limits (68.2% of the prob) by using an iterative process that moves a line down the probability axis until it intersects the curve at two points that contain 68.2% of the area between them. In short, the idea would be to adjust the horizontal black line until the red area is the desired 68.2% of the total area under the curve. Apologies for the crummy drawing.

I have absolutely no clue where to start, besides coming up with a fit to the data. I would like to make it an iterative process where I input a starting height for the horizontal line and it 'homes' in on the proper height (to within some range).

A little nudge in the right direction would be great.

Side note: I'm pretty sure this is a known technique, but I don't know what it's called. If anyone knows the answer that would be a great start.

Answer 1

[Star Strider]

I have some archive code that can help you get started.

If you want to fit what is essentially a histogram to the normal distribution, here is one way:

% Generate x- and y- data:
x = linspace(0,10);
% Normal distribution — b(1) = mean,  b(2) = std
N = @(b,x)  exp(-((x-b(1)).^2)./(2*b(2)^2)) ./ (b(2)*sqrt(2*pi));
b = [5; 1];
d1 = N(b,x);                            % Generate smooth data
d2 = d1 + 0.05*(rand(size(d1))-.5);     % Add noise

% Use ‘fminsearch’ to do regression:
y = d2;
OLS = @(b) sum((N(b,x) - y).^2);        % Ordinary Least Squares cost function
opts4 = optimset('MaxFunEvals',50000, 'MaxIter',10000);
B = fminsearch(OLS, rand(2,1))
Nfit = N(B,x);

figure(1)
plot(x, d1, '-b')
hold on
plot(x, d2, '.g')
plot(x, Nfit, '-r')
hold off
grid

I’ll let you experiment with it. Change the initial estimate for the mean to fit other peaks in your data.

With regard to the patch colouring, another bit of archived code will provide an example:

alpha = 0.05;           % significance level
mu = 82.9;              % mean
sigma = 8.698;          % std
x = linspace(mu-5*sigma, mu+5*sigma, 500);
cutoff1 = norminv(alpha/2, mu, sigma);          % Lower 95% CI is p = 0.025
cutoff2 = norminv(1-alpha/2, mu, sigma);        % Upper 95% CI is p = 0.975
y = normpdf(x, mu, sigma);—
xci = [linspace(mu-5*sigma, cutoff1); linspace(cutoff2, mu+5*sigma)];
yci = normpdf(xci, mu, sigma);

figure(1)
plot(x, y, '-b', 'LineWidth', 1.5)
patch(x, y, [0.5 0.5 0.5])
patch([xci(1,:) cutoff1], [yci(1,:) 0], [1 1 1])
patch([cutoff2 xci(2,:)], [0 yci(2,:)], [1 1 1])

This shades the 95% confidence intervals, but it easily adapted to do what you want.

[Julien]

Thanks, this helps quite a bit.

[Star Strider]

My pleasure!