| Re: Trouble with MPFITFUN [message #79944 is a reply to message #79942] |
Wed, 11 April 2012 11:08  |
Helder Marchetto
Messages: 520
Registered: November 2011
|
Senior Member |
|
|
On Wednesday, April 11, 2012 6:34:23 PM UTC+2, Helder wrote:
> Hi,
> I've been spending a bit too much time on this and I am wondering what is going wrong here.
> I'm trying to fit using a step function broadened by a Gaussian.
> The fitting function is:
>
> FUNCTION GaussStep, X, P
> ;Calculate the broadening of a step function with:
> ;P[0] = step position
> ;P[1] = left value
> ;P[2] = right value
> ;P[3] = step width
> PRINT, P
> P[0] = (P[0] > MIN(X)) < MAX(X)
> Y = DBLARR(N_ELEMENTS(X))
> LowIndeces = WHERE(X LT P[0], CountLow, COMPLEMENT = HighIndeces, NCOMPLEMENT=CountHigh)
> IF CountLow GT 0 THEN Y[LowIndeces] = P[1]
> IF CountHigh GT 0 THEN Y[HighIndeces] = P[2]
> Sigma=P[3]
> nPts=10*Sigma+1.0
> kernel=DINDGEN(nPts)-(nPts-1)/2.0
> kernel=EXP(-kernel^2/(2.*sigma^2))
> kernel/=TOTAL(kernel,/DOUBLE)
> yconvol = CONVOL(Y,kernel,/EDGE_TRUNCATE)
> RETURN, yconvol
> END
>
> To test MPFITFUN I use the following code:
> PRO TestFit
> xData = DINDGEN(201)
> yData = DBLARR(201)+RANDOMU(SEED,201,/DOUBLE)*0.2-0.1
> yData[150:200] += 1.0D
> StParam = [148D,MIN(yData),MAX(yData),3D]
> DataErr = DBLARR(N_ELEMENTS(xData))+0.2D
> Results = MPFITFUN('GaussStep', xData,yData, DataErr, StParam, STATUS=status, /quiet)
> PLOT, xData, yData
> OPLOT, xData, GaussStep(xData,Results), COLOR = 255L
> PRINT, 'Final Parameters = ', Results
> PRINT, 'Start Parameters = ', StParam
> END
>
> The output shows all the calls of the fitting function. And I find that at the end there is always NO change in the first parameter. Here is an example of the output:
>
> 148.00000 -0.099990073 1.0994661 3.0000000
> 148.00000 -0.099990073 1.0994661 3.0000000
> 148.00000 -0.099990071 1.0994661 3.0000000
> 148.00000 -0.099990073 1.0994661 3.0000000
> 148.00000 -0.099990073 1.0994661 3.0000000
> 148.00000 0.0073445709 1.0082363 2.3488363
> 148.00000 0.0073445709 1.0082363 2.3488363
> 148.00000 0.0073445710 1.0082363 2.3488363
> ...
> 148.00000 -0.0039705287 0.99188729 2.0999998
> 148.00000 -0.0039705257 0.99188729 2.1000000
> 148.00000 -0.0039705254 0.99188729 2.1000000
> 148.00000 -0.0039705254 0.99188729 2.1000000
> Final Parameters = 148.00000 -0.0039705254 0.99188729 2.1000000
> Start Parameters = 148.00000 -0.095071379 1.0978406 3.0000000
>
> Throughout all the fitting procedure the first parameter has never been changed.
>
> Am I doing something terribly wrong? I generally have no estimates for the errors in the data, therefore I used 0.1. In the example data this is easy to calculate, but the fitting has to be applied to the most different data sets.
>
> I also tried playing with the XTOL parameter without any success.
>
> Any tips are appreciated.
>
> Many thanks,
> Helder
>
> PS: I tried lots of different initial conditions, I tried using "parinfo.fixed" to block the other parameters, ... but at the end I never get any change in P[0]... sigh..
>
> PSS: The function GaussStep is working fine... I can replot the data in the correct way by moving the parameters by hand.
Ok, just had some dinner and the head is a bit clearer now. I think everything is working fine. The only problem is that the square residuals don't change a lot when changing the position of the step. It's like finding the minimum on a flat 100 square meter surface with a little 1 square cm hole hidden somewhere...
I think this is it, if anybody has a suggestion on how to better estimate the step, that would be appreciated... I'll go for dessert and maybe get a good idea on how to do that.
Cheers,
Helder
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
