Re: Re: Gradient of a List
- To: mathgroup at smc.vnet.net
- Subject: [mg82647] Re: [mg82614] Re: Gradient of a List
- From: DrMajorBob <drmajorbob at bigfoot.com>
- Date: Sat, 27 Oct 2007 06:00:27 -0400 (EDT)
- References: <10340059.1193232789565.JavaMail.root@m35> <ffpps9$l54$1@smc.vnet.net> <10446382.1193425327958.JavaMail.root@m35>
- Reply-to: drmajorbob at bigfoot.com
Scott,
I hadn't noticed I was perturbing the x-values and using UNperturbed x to
calculate y. No wonder the result looked so awful!
Your Gaussian smoothing works amazingly well in this situation.
I think I'd write it this way, however:
ClearAll[pdf, g, x, data]
pdf[xIdeal_, sig_] = Exp[-((# - xIdeal)/sig)^2/2]/(Sqrt[2 Pi] sig) &;
g[data_, sigma_] := g[data, sigma] = Function[{x},
Evaluate@Block[{x0, y0},
{x0, y0} = Transpose@data;
y0.Normalize[pdf[x, sigma] /@ x0, Total]]]
data = Table[{x + RandomReal[], Sin@x + 0.1 RandomReal[]}, {x, 0, Pi,
0.1}];
Plot[g[data, 0.5]'[x], {x, min, max}, PlotRange -> All]
Notice the "just-in-time" definition of g.
Bobby
On Fri, 26 Oct 2007 04:21:23 -0500, Scott Hemphill
<hemphill at hemphills.net> wrote:
> DrMajorBob <drmajorbob at bigfoot.com> writes:
>
>> data = Table[{x + RandomReal[], Sin@x + 0.1 RandomReal[]}, {x, 0, Pi,
>> 0.1}];
>> f = Interpolation[data, InterpolationOrder -> 3];
>> {min, max} = data[[Ordering[data][[{1, -1}]], 1]];
>> Quiet@Plot[f'[x], {x, min, max}, PlotRange -> All]
>>
>> I use Quiet because Plot sometimes samples outside the data range and
>> throws the InterpolatingFunction::dmval message.
>>
>> Notice, however, the result isn't even close to Cos[x], and it changes
>> quite a bit if you change the InterpolationOrder.
>
> Of course, these problems are because of the noise in both the x and y
> data values. Since Interpolation insists on passing exactly through
> the points given, the interpolating function has to wiggle around a
> lot to fit all the noise. The OP may not have any noise in his
> independent variables (x,y) and may have little or none in his
> function values.
>
> Still, yours is an interesting problem. One way of handling it would
> be to interpolate via weighted averages. For example, you could
> assign a Gaussian weight to all the function values based on how close
> the x value is to the x coordinates of the data:
>
> (* Gaussian centered at x0, with standard deviation sig *)
>
> pdf[x_,x0_,sig_] := 1/(Sqrt[2Pi]sig) Exp[-(x-x0)^2/(2sig^2)];
>
> (* Gaussian weighted average of data, using sig = 0.5 *)
> (* try using other values for sig *)
>
> g[x_] = Block[{x0,y0,w},
> x0 = data[[All,1]]; (* x-coordinates of data *)
> y0 = data[[All,2]]; (* y-coordinates of data *)
> w = pdf[x,#,0.5]& /@ x0; (* weight the x-coordinates *)
> w /= Plus @@ w; (* Normalize the weights *)
> w . y0 (* Return interpolated function value *)
> ];
>
> Now you have a continuous function g[x], and you can plot it as well
> as g'[x]. (Of course it is inefficient, since it recalculates all the
> weights every time you call it. You could enter "foo[x_]=g[x];" and
> then "foo[x]" wouldn't have that problem.)
>
> One nice feature of Gaussian weights happens if you have a set of
> equally spaced data points. If they extend infinitely in both
> positive and negative directions, or *equivalently* the function
> values are zero beyond the region of interest, then you can omit the
> normalization step. (One example is in image processing, where the
> region beyond the boundaries of the image may be assumed to be black.)
> Then the Gaussian weighting is equivalent to convolution(1) with a
> Gaussian kernel. This convolution has some nice properties. For
> example, it is infinitely differentiable, because the Gaussian is.
> Also, you can express its derivatives in the same form as the
> convolution itself, i.e. the convolution of a Gaussian with a set of
> data points.
>
> The OP might be able to use the two-dimensional version of the the
> Gaussian weighted interpolation above, but Mathematica's built-in
> polynomial interpolation might work perfectly well.
>
> (1) when I speak of convolution with a "data point" (x0,y0), I really
> mean convolving with the function y0*DiracDelta[x0]. The result
> is a y0 times a Gaussian centered at x0. Convolution with a
> collection of data points gives the sum of all the Gaussians.
>
>> On Wed, 24 Oct 2007 03:34:28 -0500, olalla <operez009 at ikasle.ehu.es>
>> wrote:
>>
>> > Hi everybody,
>> >
>> > Does anybody know how can I get the "gradient" of a list of points?
>> >
>> > My real problem is:
>> >
>> > I have a scalar field previously obtained numerically that for a
>> > given point (xi,yi) takes a value f(xi,yi). What I want to do is an
>> > estimation of the gradient of this scalar field BUT I haven't got any
>> > analytical function that expresses my field so I can't use the Grad
>> > function.
>> >
>> > How can I solve this using Mathematica?
>> >
>> > Thanks in advance
>> >
>> > Olalla, Bilbao UPV/EHU
>> >
>> >
>>
>>
>>
>> --
>>
>> DrMajorBob at bigfoot.com
>>
>
--
DrMajorBob at bigfoot.com