j.coenia@gmail.com wrote:
> The user function has help statements at its beginning and end now,
> and I tried the other suggestions, but I'm still stumped. I'm doing
> something foolish.
>
> Below is the problem code. Most of it is just the 353-element data
> array. Without the data, it's only two short functions, 20 lines of
> code or so total. I hope it's not bad netiquette to post so many
> lines. I can't figure out why these data crash the program. I can
> fit most other data acquired in this manner -- I have fit hundreds of
> thousands of these curves with this method, but the error crops up
> once every 50 cases or so (usually for data that don't fit the model
> very well, but it still shouldn't crash?).
>
> The data are being fit to a gamma variate function. Two things I'm
> not sure about: (1) I'm passing t0 as a fitting parameter, and I'm not
> sure if this is acceptable. (2) Worse, I don't have a good idea of
> how to estimate error, so I'm just using the standard deviation of the
> baseline (pre-arrival) curve, which is probably wrong. Suggestions?
>
> The code crashes after just two MPFIT iterations. Just type
> 'MPFIT_problem' to run.
>
> Thanks for any help.
>
> ------------------------------------
>
> function gamma_variate, x, p
>
> print, 'p start...'
> help, p
>
> baseline = p[0]
> t0 = p[1]
> tmax = p[2] ; time of peak signal
> smax = p[3] ; signal max in mVolts
> alpha = p[4]
>
> ; shift time, apply gamma variate only to data after t0
>
> wh = Where(x GE t0, ct)
I think that you may have a problem here when ct EQ 0 and wh EQ -1.
(that is, all x<t0).
Ciao,
Paolo
> t = (x - t0)[wh]
> tmax = tmax - t0
>
> ; gamma variate function
>
> s = baseline + (smax) * (tmax^(-alpha)) * exp(alpha) * (t^alpha) * exp
> ((-alpha) * t / tmax)
>
> ; add the pre-t0 data (baseline)
>
> n = n_elements(x)
> pre = make_array(n - ct, Value = baseline)
> s = [pre, s]
>
> print, 'p end...'
> help, p
>
> return, s
>
> end;-------------------------------------------------------- --
>
>
> pro mpfit_problem
>
> !except = 2
>
> data = [ $
> 64.1682 , $
> 66.3804 , $
> 73.4243 , $
> 64.1682 , $
> 64.1682 , $
> 55.9551 , $
> 47.0046 , $
> 52.2084 , $
> 48.6855 , $
> 52.2084 , $
> 57.9162 , $
> 75.9137 , $
> 64.1682 , $
> 75.9137 , $
> 86.6244 , $
> 89.4994 , $
> 89.4994 , $
> 101.845 , $
> 105.153 , $
> 89.4994 , $
> 95.4993 , $
> 92.457 , $
> 105.153 , $
> 89.4994 , $
> 101.845 , $
> 92.457 , $
> 86.6244 , $
> 83.8301 , $
> 89.4994 , $
> 92.457 , $
> 92.457 , $
> 95.4993 , $
> 95.4993 , $
> 95.4993 , $
> 89.4994 , $
> 89.4994 , $
> 75.9137 , $
> 71.0068 , $
> 78.4766 , $
> 57.9162 , $
> 64.1682 , $
> 59.9377 , $
> 52.2084 , $
> 45.3754 , $
> 39.3496 , $
> 39.3496 , $
> 48.6855 , $
> 62.0212 , $
> 73.4243 , $
> 50.4197 , $
> 73.4243 , $
> 66.3804 , $
> 66.3804 , $
> 68.6594 , $
> 68.6594 , $
> 86.6244 , $
> 57.9162 , $
> 89.4994 , $
> 73.4243 , $
> 73.4243 , $
> 78.4766 , $
> 75.9137 , $
> 83.8301 , $
> 95.4993 , $
> 86.6244 , $
> 73.4243 , $
> 81.1148 , $
> 86.6244 , $
> 89.4994 , $
> 40.785 , $
> 54.053 , $
> 68.6594 , $
> 78.4766 , $
> 73.4243 , $
> 95.4993 , $
> 83.8301 , $
> 81.1148 , $
> 95.4993 , $
> 83.8301 , $
> 75.9137 , $
> 57.9162 , $
> 68.6594 , $
> 83.8301 , $
> 78.4766 , $
> 75.9137 , $
> 68.6594 , $
> 78.4766 , $
> 75.9137 , $
> 83.8301 , $
> 75.9137 , $
> 64.1682 , $
> 71.0068 , $
> 78.4766 , $
> 73.4243 , $
> 64.1682 , $
> 92.457 , $
> 68.6594 , $
> 73.4243 , $
> 75.9137 , $
> 62.0212 , $
> 73.4243 , $
> 55.9551 , $
> 57.9162 , $
> 62.0212 , $
> 73.4243 , $
> 75.9137 , $
> 75.9137 , $
> 92.457 , $
> 83.8301 , $
> 112.047 , $
> 98.6279 , $
> 123.118 , $
> 92.457 , $
> 105.153 , $
> 92.457 , $
> 75.9137 , $
> 86.6244 , $
> 92.457 , $
> 73.4243 , $
> 78.4766 , $
> 66.3804 , $
> 75.9137 , $
> 119.327 , $
> 127.011 , $
> 105.153 , $
> 73.4243 , $
> 54.053 , $
> 73.4243 , $
> 75.9137 , $
> 66.3804 , $
> 86.6244 , $
> 68.6594 , $
> 78.4766 , $
> 75.9137 , $
> 78.4766 , $
> 81.1148 , $
> 127.011 , $
> 152.651 , $
> 95.4993 , $
> 83.8301 , $
> 119.327 , $
> 50.4197 , $
> 68.6594 , $
> 95.4993 , $
> 75.9137 , $
> 131.009 , $
> 54.053 , $
> 108.553 , $
> 62.0212 , $
> 83.8301 , $
> 57.9162 , $
> 152.651 , $
> 119.327 , $
> 119.327 , $
> 54.053 , $
> 73.4243 , $
> 86.6244 , $
> 89.4994 , $
> 105.153 , $
> 73.4243 , $
> 78.4766 , $
> 112.047 , $
> 50.4197 , $
> 89.4994 , $
> 92.457 , $
> 30.5028 , $
> 55.9551 , $
> 119.327 , $
> 131.009 , $
> 68.6594 , $
> 71.0068 , $
> 81.1148 , $
> 112.047 , $
> 73.4243 , $
> 64.1682 , $
> 52.2084 , $
> 123.118 , $
> 127.011 , $
> 75.9137 , $
> 112.047 , $
> 101.845 , $
> 95.4993 , $
> 105.153 , $
> 55.9551 , $
> 135.114 , $
> 83.8301 , $
> 95.4993 , $
> 54.053 , $
> 36.6134 , $
> 42.2669 , $
> 71.0068 , $
> 68.6594 , $
> 95.4993 , $
> 115.638 , $
> 57.9162 , $
> 148.095 , $
> 48.6855 , $
> 75.9137 , $
> 95.4993 , $
> 127.011 , $
> 187.992 , $
> 112.047 , $
> 172.084 , $
> 71.0068 , $
> 92.457 , $
> 135.114 , $
> 148.095 , $
> 36.6134 , $
> 89.4994 , $
> 105.153 , $
> 68.6594 , $
> 119.327 , $
> 27.2923 , $
> 55.9551 , $
> 95.4993 , $
> 83.8301 , $
> 152.651 , $
> 43.7966 , $
> 62.0212 , $
> 112.047 , $
> 123.118 , $
> 135.114 , $
> 105.153 , $
> 32.8277 , $
> 62.0212 , $
> 39.3496 , $
> 47.0046 , $
> 92.457 , $
> 108.553 , $
> 86.6244 , $
> 152.651 , $
> 55.9551 , $
> 119.327 , $
> 52.2084 , $
> 75.9137 , $
> 71.0068 , $
> 35.3102 , $
> 43.7966 , $
> 119.327 , $
> 112.047 , $
> 127.011 , $
> 127.011 , $
> 62.0212 , $
> 54.053 , $
> 177.256 , $
> 311.315 , $
> 123.118 , $
> 101.845 , $
> 92.457 , $
> 152.651 , $
> 127.011 , $
> 101.845 , $
> 68.6594 , $
> 162.121 , $
> 172.084 , $
> 167.04 , $
> 152.651 , $
> 108.553 , $
> 172.084 , $
> 143.655 , $
> 152.651 , $
> 162.121 , $
> 279.383 , $
> 162.121 , $
> 123.118 , $
> 34.0487 , $
> 62.0212 , $
> 83.8301 , $
> 177.256 , $
> 172.084 , $
> 157.326 , $
> 73.4243 , $
> 135.114 , $
> 92.457 , $
> 59.9377 , $
> 131.009 , $
> 42.2669 , $
> 71.0068 , $
> 78.4766 , $
> 135.114 , $
> 78.4766 , $
> 143.655 , $
> 112.047 , $
> 78.4766 , $
> 105.153 , $
> 66.3804 , $
> 73.4243 , $
> 89.4994 , $
> 143.655 , $
> 75.9137 , $
> 152.651 , $
> 89.4994 , $
> 98.6279 , $
> 95.4993 , $
> 66.3804 , $
> 24.3893 , $
> 86.6244 , $
> 57.9162 , $
> 119.327 , $
> 135.114 , $
> 75.9137 , $
> 40.785 , $
> 112.047 , $
> 47.0046 , $
> 89.4994 , $
> 92.457 , $
> 95.4993 , $
> 89.4994 , $
> 50.4197 , $
> 81.1148 , $
> 57.9162 , $
> 45.3754 , $
> 66.3804 , $
> 73.4243 , $
> 59.9377 , $
> 73.4243 , $
> 105.153 , $
> 78.4766 , $
> 64.1682 , $
> 54.053 , $
> 62.0212 , $
> 39.3496 , $
> 47.0046 , $
> 89.4994 , $
> 28.327 , $
> 81.1148 , $
> 43.7966 , $
> 68.6594 , $
> 48.6855 , $
> 54.053 , $
> 40.785 , $
> 35.3102 , $
> 42.2669 , $
> 75.9137 , $
> 64.1682 , $
> 64.1682 , $
> 39.3496 , $
> 95.4993 , $
> 78.4766 , $
> 101.845 , $
> 64.1682 , $
> 68.6594 , $
> 48.6855 , $
> 59.9377 , $
> 66.3804 , $
> 83.8301 , $
> 40.785 , $
> 39.3496 , $
> 78.4766 , $
> 127.011 , $
> 172.084 , $
> 52.2084 , $
> 119.327 , $
> 37.9595 , $
> 66.3804 , $
> 108.553 , $
> 86.6244 , $
> 83.8301 $
> ]
>
> ; calculate times @ 30Hz
>
> n = n_elements(data)
> t = findgen(n)/30.0
>
> ; define and contrain parameters
>
> parinfo = replicate({value:0.D, fixed:0, limited:[0,0], $
> limits:[0.D,0], mpmaxstep: 0.0D}, 5)
>
> ; baseline
> parinfo(0).fixed = 0b
> parinfo(0).limited = [1,1]
> parinfo(0).limits = [0, 3675.36]
> parinfo(0).value = 84.551201
>
> ; t0
> t0 = 4.0333333
> parinfo(1).fixed = 0b
> parinfo(1).limited(0) = 1
> parinfo(1).limits(0) = 0.D
> parinfo(1).value = 4.0333333
>
> ; tmax
> parinfo(2).fixed = 0b
> parinfo(2).limited = [1,1]
> parinfo(2).limits = [t0, n/30.0 - 0.5] ; 0.5 secs before end
> parinfo(2).value = 8.5333338
>
> ; max
> parinfo(3).limited = [1,1]
> parinfo(3).limits = [0, 3675.36]
> parinfo(3).value = 127.96773
>
> ; alpha
> parinfo(4).fixed = 0b
> parinfo(4).limited = [1,1]
> parinfo(4).limits = [0.5, 12]
> parinfo(4).value = 1.0
>
> ; calculate std deviation of signal up to t0
> ; probably not the best way to estimate error...
>
> err = stddev(data[0:t0*30]) > 1.0
>
> ; plot the signal
>
> window, /free
> plot, t, data, $
> Title = 'Signal', $
> XTitle = 'Time (s)', $
> YTitle = 'mVolts'
>
> ; fit
>
> parms = mpfitfun('gamma_variate', t, data, replicate(err, n), $
> YFit = fit, $
> Bestnorm = chisq, $
> /NaN, $
> Parinfo = parinfo, $
> Quiet = 0 $
> )
>
> ; plot fit
>
> oplot, t, gamma_variate(t, parms), Color = FSC_Color('red')
>
> end;-------------------------------------------------------- --
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
