MATLAB Answers

K
0

Fit and Plot Gaussian Function

Asked by K
on 24 Jul 2013
I'm trying to fit a Gaussian function to my data.
I use
pd=fitdist(y,'Normal');
x_values = -250:1:250;
p = pdf(pd,x_values);
plot(x_values,p)
to fit and plot the function.
It looks to be the right shape, however, the function itself is very small (the max only coming to about 4*10^-3). I know that a normal function dictates that the integral go to 1, but is there any way to keep the shape, just make it bigger so that it can plot on top of my data (X range -200, 200 Y range -250, 250)? Or is there another kind of function that will have the same shape but forgo the integration to one part?
Thanks in advance.

  1 Comment

dpb
on 24 Jul 2013
Might try using
dfittool
instead. I've not played w/ the fitdist() thingie yet; I'm somewhat puzzled regarding how it functions on a single vector--and the doc is certainly unclear on what x is really expected to be at first blush, anyway. Is it an order vector--I'm not sure otomh.
Certainly unless the estimate of sigma is quite large the standardized z of (+/-250-mu)/sigma is going to be quite large and the resulting pdf magnitudes very small indeed.
Presuming mu to be ~1 given the range is symmetric in x_values, then unless sigma is roughly 85 or larger, z starts at something >3. For that the distribution will be pretty broad.
Alternatively, for plotting purposes, you can scale either by the ratio of the two at a given (set of) points. If the x-vector is a cdf order vector as I guess must be then the midpoint magnitude should match that of the estimated mean z=0 point, also the 50% cdf point.

Sign in to comment.

1 Answer

Answer by Richard Brown on 24 Jul 2013
Edited by Richard Brown on 24 Jul 2013

Sure, rather than trying to fit a distribution (which is not what you want), just fit the Gaussian itself.
Generate some noisy Gaussian data:
a = 5; b = 17; c = 40;
noise_sd = .1;
n = 100;
x = linspace(-250, 250, n);
y = a*exp(-(x-b).^2/c^2) + noise_sd * randn(size(x));
Define the parametrised functional form.
f = @(x, p) p(1) * exp(-((x - p(2))/p(3)).^2);
And solve using lsqnonlin. You might need to be a little bit careful with initial conditions:
p = lsqnonlin(@(p) f(x, p) - y, [1 1 1]);
plot(x, y, 'ro', x, f(x, p))

  0 Comments

Sign in to comment.