Fit with multiple gaussians in multi-dimension
MG
what would be to best approach in writing a fit function that is the sum of an arbitrary number of gaussian functions in arbitrary dimensions (including co-variance).
I would try to work with vectors/matrices as much as possible, so this code should work with any dimension (assuming waves have the correct DimSize).
//multidim gauss with vectors and covariance matrix
//wVectorMu contains center of gaussian
//wMatrixCoVar is co-variance matrix
//wVectorVar is vector of variables, i.e. {x1, x2, x3}
MatrixOP /FREE intermediate1 = exp(-1/2 * ( (wVectorVar -wVectorMu)^t x Inv(wMatrixCoVar) x (wVectorVar -wVectorMu) ))
MatrixOP /FREE intermediate2 = intermediate1 * 1/(sqrt( powR(2*pi, numdims) * Det(wMatrixCoVar) ))
returnVal = intermediate2[0]
//wVectorMu contains center of gaussian
//wMatrixCoVar is co-variance matrix
//wVectorVar is vector of variables, i.e. {x1, x2, x3}
MatrixOP /FREE intermediate1 = exp(-1/2 * ( (wVectorVar -wVectorMu)^t x Inv(wMatrixCoVar) x (wVectorVar -wVectorMu) ))
MatrixOP /FREE intermediate2 = intermediate1 * 1/(sqrt( powR(2*pi, numdims) * Det(wMatrixCoVar) ))
returnVal = intermediate2[0]
What is still missing is how to sum up these gaussians with weighing and how to put all the parameters (wWeight, wVectorMu, wMatrixCoVar [there are # of gaussian different ones for the later two]) into a FitFunc. As far as I understand it, only one parameter wave is allowed in functions of type FitFunc.
Cheers,
MG
You could write a fit function that looks at the number of fit coefficients in the coefficient wave to compute the number of Gaussians to fit.
For the 2D case, if you are feeling adventurous, you could try the sum-of-fit-functions approach, and use the built-in Gauss2D fit function. If you do this, you need to hold the Y0 coefficient for all but one of the Gaussians.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
February 12, 2016 at 10:02 am - Permalink
I somehow got this to work for 3D & arbitrary number of gaussian functions (called "states" in the code). The FitFunc looks like this:
Function FitMultiGauss3D(w, x1, x2, x3) : FitFunc
WAVE w
Variable x1
Variable x2
Variable x3
Variable returnValue = 0
Variable state
Variable NumStates = numpnts(w)/13
Variable NumDims = 3
if (mod(numpnts(w), 13) != 0)
Print "Check length of W_coef."
Abort
endif
//CustomFitDialog/ Independent Variables 3
//CustomFitDialog/ x1
//CustomFitDialog/ x2
//CustomFitDialog/ x3
//CustomFitDialog/ Coefficients NumState*13
//CustomFitDialog/ w[0+12*state,11+12*state] = wVectorMu, wMatrixCoVar
//CustomFitDialog/ w[12*NumStates,12*NumStates+3] = wPi
Duplicate /O/R=[0,12*NumStates-1] w, b_params
Redimension /N=(3,4,NumStates) b_params
Duplicate /O /R=[12*NumStates,12*NumStates+3] w, wPi
for (state=0; state<NumStates; state+=1)
//extract means Mu and covariance matrix
Duplicate /FREE /R=[0,*][0][state] b_params, wVectorMu
Redimension /N=(-1) wVectorMu
Duplicate /FREE /R=[0,*][1, *][state] b_params, wMatrixCoVar
Redimension /N=(-1, -1) wMatrixCoVar
//Make wVectorVar
Make /O/N=3 wVectorVar={x1, x2, x3}
//multidim gauss with vectors and covariance matrix
//wVectorMu contains center of gaussian
//wMatrixCoVar is co-variance matrix
//wVectorVar is vector of variables, i.e. {x1, x2, x3}
MatrixOP /FREE intermediate1 = exp(-1/2 * ( (wVectorVar -wVectorMu)^t x Inv(wMatrixCoVar) x (wVectorVar -wVectorMu) ))
MatrixOP /FREE intermediate2 = intermediate1 * 1/(sqrt( powR(2*pi, NumDims) * Det(wMatrixCoVar) ))
returnValue += intermediate2[0] * wPi[state]
endfor
return returnValue
End
The W_coef has a fixed format how the 13 parameters per gaussian function are sorted. This helps extracting the vector of the center of the gaussian (wVectorMu) and the covariance matrix by just Redimension the appropriate range of the W_coef.
The fit is called by
FuncFitMD /N /W=2 /NTHR=0 FitMultiGauss3D, W_coef, DataWave, but it is painfully slow.February 29, 2016 at 10:11 am - Permalink