Creating an arbitrary number of displaced Gaussian

Hi, I have created a set of 3 gaussians as follows:

%add 3 gaussians together
a=100;a1=75;a2=123
b=10;b1=25;b2=40
c=1.5;c1=1.5;c2=1.5
d=1;d1=3.2;d2=2

x = 0:1:100;
gaussEqn1 =a*exp(-0.5*((x-b)/(c)).^2)+d+a1*exp(-0.5*((x-b1)/(c1)).^2)+d1+a2*exp(-0.5*((x-b2)/(c2)).^2)+d2

This gives me:

Is it possible to redefine the gaussian equation and coefficients a,b,c,d so I can add an arbitrary number of gaussians - I will need upto 20 but would like to have the ability of any number upto 20. I would like to constrain the parameters such as

a=between 1000 and 2000 - randomly
b=15, then the next b is 30, then the next is 45 so the spacing between them if fixed.
c=1.5 for all gaussians
d can vary randomly between 0 - 10.

Thanks Jason

Answer 1

[Image Analyst]

Just put your code into a loop:

c = 1.5;
numGaussians = 10;
x = linspace(0, numGaussians * 16, 300);
gaussEqn1 = zeros(1, length(x));
for k = 1 : numGaussians
  % Construct this one Gaussian
  a = 1000 + 1000 * rand(1);
  b = 15 * k;
  d = 10 * rand(1);
  fprintf('For k = %d, a = %f, b = %f, c = %.1f, d = %f\n', k, a, b, c, d);
  thisGaussian = a * exp(-0.5 * ((x-b)/c) .^ 2) + d;
  % Add into accumulator array:
  gaussEqn1 = gaussEqn1 + thisGaussian;
end
plot(x, gaussEqn1, 'b-', 'LineWidth', 2);
grid on;

Answer 2

[David J. Mack]

Hey Jason, BSXFUN is your friend (and your enemy w.r.t. readability):

%Parameters - must be scalar or 1xm row vectors.
a=[100 75 123];  % or a=1000*rand(1,m)+1000;
b=[10 25 40];    % or b=15:15:15*m;
c=[1.5 1.5 1.5]; % or c=1.5;
d=[1 3.2 2];     % or d=10*rand(1,m);
%x must be nx1 column vector.
x = (0:1:100)';

% Function handle takes arbitrary inputs, creates matrix of Gaussians with the specified parameters
gaussFun = @(x,a,b,c,d)bsxfun(@plus,bsxfun(@times,exp(-0.5*bsxfun(@rdivide,bsxfun(@minus,x, b), c).^2), a), d);

What does it do: This is exactly the same formula as yours, except that BSXFUN takes care of the non-matching dimension by expanding them accordingly. E.g. (nx1).*(1xm) -> (nxm) but not in a mathematical sense (e.g. matrix multiplication).

You would use it to create the sum-of-gaussians in the following way:

gaussEqn1 = sum(gaussFun(x,a,b,c,d),2);

The function handle can easily be adapted to a different domain, e.g. x2=(-50:50)' (although that doesn't make too much sense).

gaussEqn2 = sum(gaussFun(x2,a,b,c,d),2);

Or you can visualize each single Gaussian by:

plot(x,gaussFun(x,a,b,c,d));

With this formulation you can use arbitrarily sized parameter vectors (each with the same size) generating a column in the resulting nxm output of the function handle.

Greetings, David

[Jason]

David, thanks for your answer. Sorry I didn't accept your answer, but I.A's was simpler to follow. But I am interested to take some time to follow yours through.

[David J. Mack]

That's probably the fate of a high-level answer...