4
$\begingroup$

I am attempting to fit to the convolution of a simple function with a Gaussian. (I realize this could be done with Convolve also, but I thought I could make it go faster this way). I am able to define the function itself and the convolved version with

ExactFunction[x_, w_, h_, x0_] := h*UnitStep[x + w/2 - x0] UnitStep[w/2 - x + x0]
Gaussian[x_, δ_] := 1/(δ Sqrt[2 π]) E^(-x^2/(2 δ^2))
ConvolvedFunction[x_?NumericQ, w_?NumericQ, h_?NumericQ, x0_?NumericQ, δ_?NumericQ] := 
 NIntegrate[ ExactFunction[x - xp, w, h, x0]* Gaussian[xp, δ], {xp, -5 δ, 5 δ},  AccuracyGoal -> 4]
Show[Plot[ExactFunction[x, 50, 2, 50], {x, 0, 100}, PlotStyle -> Dashed], 
     Plot[ConvolvedFunction[x, 50, 2, 50, 5], {x, 0, 100}]]

where I have used ?NumericQ to prevent ConvolvedFunction from evaluating prematurely. I can fit to the original function easily with something like

Data = Table[Sin[x π/100], {x, 0, 100}] + RandomReal[.2, 101];
FitParameters = FindFit[Data, ExactFunction[x, w, h, x0], {{w, 50}, {h, 1}, {x0, 50}}, x, 
                        Method -> "PrincipalAxis"]
Show[ListPlot[Data],  Plot[ExactFunction[x, w, h, x0] /. FitParameters, {x, 0, 100}]]

However, if I try to fit the ConvolvedFunction even for small MaxIterations it never finishes.

ConvolvedFitParameters = FindFit[Data, ConvolvedFunction[x, w, h, x0, δ], 
            {{w, 81}, {h, .86}, {x0, 52}, {δ, 5}}, x, Method -> "PrincipalAxis", MaxIterations -> 10]
Show[ListPlot[Data], Plot[ConvolvedFunction[x, w, h, x0, δ] /. ConvolvedFitParameters, {x, 0, 100}]]

What have I done wrong? I've looked at a number of similar examples and don't see any obvious differences. Thanks for any help you can provide.

Edit: The solution provided by @belisarius works in the simple example I posted, but I think it is using the fact that the convolution of a Gaussian with a step function has a closed form. My actual function is

VSumExact[x_, w_, CSum_, dposdx_, x0_ ] :=   
CSum/w (Log[1/4 (w + 2 dposdx (x - x0))^2] -
Log[1/4 (w - 2 dposdx (x - x0))^2])*UnitStep[w/2 - dposdx (x - x0)] 
UnitStep[dposdx (x - x0) + w/2]

and when I try to run

VSumSmooth[x_, w_, CSum_, dposdx_, x0_ , \[Delta]_] :=  
Evaluate@Convolve[Gaussian[xp, \[Delta]], 
VSumExact[xp, w, CSum, dposdx, x0], xp, x]

it never finished executing.

$\endgroup$
3
$\begingroup$
cf1[x_, w_, h_, x0_, δ_] :=  Evaluate@Convolve[Gaussian[xp, δ],  ExactFunction[xp, w, h, x0], xp, x]
ConvolvedFitParameters = FindFit[Data, cf1[x, w, h, x0, δ], {{w, 50}, {h, 1}, {x0, 50}, {δ, 5}}, x]
Show[ListPlot[Data], Plot[cf1[x, w, h, x0, δ] /. ConvolvedFitParameters, {x, 0, 100}]]

Mathematica graphics

$\endgroup$
  • $\begingroup$ Welp, that works, thanks! I'll read up on Evaluate@. $\endgroup$ – ARM Jan 24 '15 at 18:32
  • $\begingroup$ Well, the elation was short lived -- in my actual code the function I am trying to convolve is a bit more complicated, and I can't even plot it after a few min, let alone fit to it. I'll try optimizing it or writing the NIntegrate version in a similar form... $\endgroup$ – ARM Jan 24 '15 at 18:43
3
$\begingroup$

Here's how I eventually got this working with my actual function. I think what I originally posted was correct (or close?), I just needed to set AccuracyGoal low enough. Hopefully this helps anyone in a similar situation.

Setting up the convolved function:

BeamProfile[x_, \[Delta]_] := 
 1/(\[Delta] Sqrt[2 \[Pi]]) E^(-x^2/(2 \[Delta]^2))
VSumExact[x_, w_, CSum_, dposdx_, x0_ ] := 
 CSum/w (Log[1/4 (w + 2 dposdx (x - x0))^2] - 
 Log[1/4 (w - 2 dposdx (x - x0))^2])*
 UnitStep[w/2 - dposdx (x - x0)] UnitStep[dposdx (x - x0) + w/2]
VSumSmooth[x_?NumericQ, w_?NumericQ, CSum_?NumericQ, dposdx_?NumericQ,
 x0_ ?NumericQ, \[Delta]_?NumericQ] := 
NIntegrate[VSumExact[x - xp, w, CSum, dposdx, x0]*
 BeamProfile[xp, \[Delta]], {xp, -5 \[Delta], 5 \[Delta]}, 
 AccuracyGoal -> 4]
Show[Plot[VSumExact[x, 50*10^-6, 5*10^-10, 2*10^-6, 25], {x, 0, 50}, 
 PlotRange -> All, PlotStyle -> Dashed],
 Plot[VSumSmooth[x, 50*10^-6, 5*10^-10, 2*10^-6, 25, 1], {x, 0, 50}, 
 PlotRange -> All, PlotPoints -> 50]]

enter image description here

Making initial guesses and running the fit:

Data = {-2.19919`*^-6, -7.23919`*^-7, -4.74143`*^-7, -3.80621`*^-7, \
 -1.00522`*^-6, -4.3233`*^-7, -1.15395`*^-6, -1.27838`*^-6, \
 -5.20207`*^-7, 7.19011`*^-7, 9.42227`*^-6, 0.0000329935`, 
 0.0000454738`, 0.0000394872`, 0.0000327209`, 0.0000260497`, 
 0.0000233109`, 0.0000195878`, 0.000016143`, 0.0000123143`, 
 9.57289`*^-6, 5.64377`*^-6, 2.82322`*^-6, 
 1.60046`*^-6, -1.14659`*^-6, -4.04984`*^-6, -7.78951`*^-6, \
 -0.000010962`, -0.0000138128`, -0.000015891`, -0.0000217059`, \
 -0.0000234981`, -0.0000276923`, -0.0000331199`, -0.0000422404`, \
 -0.000044107`, -0.0000282739`, -6.27234`*^-6, -1.18261`*^-6, \
 -1.71811`*^-6, -1.54246`*^-6, -2.24305`*^-6, -1.64444`*^-6, \
 -6.59291`*^-7, -2.06042`*^-6};
width = 50*10^-6; x0guess = 25; dposdxguess = 
 2*10^-6; \[Delta]guess = 1; CSumguess = -5*10^-10;
Show[ListPlot[Data, Mesh -> All, 
 AxesLabel -> {"Position (a.u.)", "Signal (V)"}], 
Plot[VSumSmooth[x, width, CSumguess, dposdxguess, 
 x0guess, \[Delta]guess], {x, 0, 50}, PlotPoints -> 50], 
 PlotRange -> All]
SmoothFit = 
 FindFit[Data, 
 VSumSmooth[x, width, CSum, dposdx, 
 x0, \[Delta]], {{CSum, CSumguess}, {dposdx, dposdxguess}, {x0, 
 x0guess}, {\[Delta], \[Delta]guess}}, x, AccuracyGoal -> 4]
Show[Plot[
 VSumSmooth[x, width, CSum, dposdx, x0, \[Delta]] /. SmoothFit, {x, 
 0, 50}], 
ListPlot[Data, Mesh -> All, 
 AxesLabel -> {"Position (a.u.)", "Signal (V)"}]]

enter image description here

{CSum->-4.06001*10^-10,dposdx->2.00797*10^-6,x0->24.4954,\[Delta]->1.08409}

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.