1

I have a function that is a summation of several Gaussians. Working with a 1D Gaussian here, there are 3 variables for each Gaussian: A, mx, and sigma:

$A \exp \left ( - \frac{\left ( x - mx \right )^{2}}{2 \times sigma^{2}} \right )$

A*Exp[-((x - mx)^2/(2 sigma^2))]

The number of Gaussians in the final function will be vary each time the function is called, so my question is: What is the best way to define a function in Mathematica that can handle this variation, rather than hard-coding each Gaussian?

I was thinking along the lines of providing a list of {A,mx,sigma} to the function, so that if I want one Gaussian, I provide:

f[{{A,mx,sigma}}]

And if I want two Gaussians, I provide

f[{{A,mx,sigma},{A2,mx2,sigma2}}]

which would give:

A*Exp[-((x-mx)^2/(2sigma^2))] + A2*Exp[-((x-mx2)^2/(2sigma2^2))]

and so on.

But I'm not at all sure how to design the function f[] to do this efficiently (for example, can it be done without a For[] loop? Can it be compiled in future if necessary?).

Any help much appreciated - I did several searches on here and couldn't find anything, but I realise that might be because I'm not sure how to define my problem succinctly, so apologies if it has been asked before and I've missed it.

m_goldberg
  • 107,779
  • 16
  • 103
  • 257
dr.blochwave
  • 8,768
  • 3
  • 42
  • 76
  • It is not clear to me that this question is statistically valid. You have a random variable $X_i \sim N(\mu_i, \sigma_i^2)$ with pdf $f_i(x_i)$. Why are you adding the pdfs $f_i(x_i)$? Are you aware that your density has to integrate to 1 in order to be well-defined? Are you intending to find the pdf of $X_1 + X_2 + ...$? Or something else? It is usually better to first express your problem as a math stats problem, before diving into the mma code. – wolfies Jun 23 '14 at 15:44
  • It's for some Gaussian peak-fitting. – dr.blochwave Jun 23 '14 at 15:46
  • Or rather, some peak-fitting that happen to look a lot like Gaussians. – dr.blochwave Jun 23 '14 at 15:47
  • 1
    For fitting Gaussians, you might be interested in the answers here. – KennyColnago Jun 23 '14 at 16:18
  • Yep, thanks Kenny - I'd seen that before (and in fact looking again it does actually include an answer to this question...oops!) and am using some of the ideas from it. – dr.blochwave Jun 23 '14 at 21:45

3 Answers3

4
f[data_] := Total[#1*Exp[-((x - #2)^2/(2 #3^2))] & @@@ data];
f[{{A, mx, sigma}}]
f[{{A, mx, sigma}, {A2, mx2, sigma2}}]

Use Function, Apply and Total.

enter image description here

Apple
  • 3,663
  • 16
  • 25
4

Since all of the component functions are Listable

f[{a_, m_, s_}, x_] := Total[a*Exp[-(x - m)^2/(2*s^2)]]

n = 5;

amp = Array[a, n];

mean = Array[m, n];

sigma = Array[s, n];

f[{amp, mean, sigma}, x]

a[1]/E^((x - m[1])^2/ (2*s[1]^2)) + a[2]/E^((x - m[2])^2/ (2*s[2]^2)) + a[3]/E^((x - m[3])^2/ (2*s[3]^2)) + a[4]/E^((x - m[4])^2/ (2*s[4]^2)) + a[5]/E^((x - m[5])^2/ (2*s[5]^2))

Bob Hanlon
  • 157,611
  • 7
  • 77
  • 198
1

The solution by Apple was very helpful. However, I found it did not work for me (I had a similar situation). No function of x is created in that example. I found the following would work:

f[data_] := Total[#1*Exp[-((x - #2)^2/(2 #3^2))] & @@@ data];
(* a specific example might be *)
g[x_] := Evaluate[f[{{x, 1.0, 0.5, 0.5}, {x, 2.0, 0.7, 1.5}}]]

g now acts as a function of x, but of course the parameters needed to be specified before numerical values are obtained. This is done above.

Anthony Mannucci
  • 463
  • 1
  • 4
  • 7