MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Graphics export to BMP and JPEG
  • Next by Date: Re: Tooltip does not what I expect it to do
  • Previous by thread: Re: Gradient of a List
  • Next by thread: Re: Re: Re: Gradient of a List