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

Re: ListIntegrate3D ?

• To: mathgroup at smc.vnet.net
• Subject: [mg76832] Re: ListIntegrate3D ?
• From: Peter Pein <petsie at dordos.net>
• Date: Mon, 28 May 2007 01:01:45 -0400 (EDT)
• References: <f38r74$i1f$1@smc.vnet.net>

nottrevorgraham at gmail.com schrieb:
> Hi,
> I have a list of x,y,z data:
> x1, y1, z1
> x2, y2, z2
> ....
> xN, yN, zN
> 
> That defines a surface z=f(x,y).  I would like to be able to integrate
> to compute the area under this surface.
> 
> To do this, I assume I need first need to interpolate the surface,
> then integrate the interpolated approximation (which will be an
> approximation to f(x,y)).  In 2D this seems straightforward, using
> NumericalMath`ListIntegrate`.  Could anyone suggest how I could do
> this in 3D?
> 
> The data is two parameters (x,y) and z is the (log)-likelihood of a
> model given the parameters.  I'm trying to integrate the likelihoods
> over (x,y) to compute a posterior distribution.
> 
> Many thanks for your help.
> 
> Trevor
> 
> 

Hello Trevor,

in the worst case your (x,y)-points do not lie on a rectangular grid. In such
a case, I use the IMTEK Mathematica Supplement; you can find it at
http://www.imtek.uni-freiburg.de/simulation/mathematica/IMSweb/ .

After installation get imsInterpolation:

In[1]:=
<< "Imtek`Interpolation`"
<< "NumericalMath`NIntegrateInterpolatingFunct`"

Here we get a disturbed array of points:

In[3]:= SeedRandom[1];
data = Table[({#1, #2, (1 + (Random[] - 1/2)/50)*
   Exp[-((1/3 - #1)^2 + (#2 - 1/2)^2)/2]} & )[4*Random[] - 2, 4*Random[] - 2],
 {70}];

now build the interpolating function:

In[5]:=
f = imsUnstructuredInterpolation[data, imsSpline[x, xi, 3]];

Say, you want to integrate over a disc, centered at place of the maximum of f
and with a radius such that f take apprimately the value 1/2 on the boundary
of the disk (because the outcome is "nice").

Get the center and the radius:

In[6]:=
{x0, y0} = {x, y} /. Last[NMaximize[{f[x, y], 0 < x < 1 > y > 0}, {x, y}]]
getr0[(t_)?NumericQ] := r /.
   FindRoot[f @@ ({x0, y0} + r*{Cos[t], Sin[t]}) == 1/2, {r, 1}];
r0 = (1/(2*Pi))*NIntegrate[getr0[t], {t, -Pi, Pi}]

Out[6]=
{0.33819365081907227, 0.5860641465257586}
Out[8]=
1.1739273460216542

just to visualize the current situation:

In[9]:=Export["viz.png",
DensityPlot[f[x, y], {x, -2, 2}, {y, -2, 2}, PlotPoints -> 65, Mesh -> False,
  AspectRatio -> Automatic, Background -> Black, Epilog -> Flatten[
    {Red, Circle[{x0, y0}, r0], AbsolutePointSize[5],
     Apply[Function[{x, y},{(#1*(1 - f[x, y])^2 & ) /@ RGBColor[1, 0.5, 0.25],
        Point[{x, y}]}], Most /@ data, {1}], Green, Point[{x0, y0}],
     Line[{{x0, y0}, {x0, y0} + r0/Sqrt[2]}]}], ImageSize -> 500]]

[ graphics of course ommitted; you can find it at
http://freenet-homepage.de/Peter_Berlin/Mathe/viz.png ]
Out[9]= "viz.png"

And here we go:

In[10]:=
NIntegrateInterpolatingFunction[f[x, y], {x, x0 - r0, x0 + r0},
  {y, y0 - Sqrt[r0^2 - (x - x0)^2], y0 + Sqrt[r0^2 - (x - x0)^2]},
  Method -> DoubleExponential]
Out[10]=
3.079901416226034

The input has an error-range of 2% and so an error of
In[11]:=
100*(Abs[% - Pi]/Pi)
Out[11]=
1.963693074379535

percent is quite good, isn't it?

(just to verify, that Pi should be the outcome:
In[12]:=
Block[{r0 = Sqrt[2*Log[2]], x0 = 1/3, y0 = 1/2,
 f = Exp[-((1/3 - #1)^2 + (#2 - 1/2)^2)/2] & },
  NIntegrate[f[x, y], {x, x0 - r0, x0 + r0},
   {y, y0 - Sqrt[r0^2 - (x - x0)^2], y0 + Sqrt[r0^2 - (x - x0)^2]},
   Method -> DoubleExponential]
] - Pi
Out[12]=
0.
).

Regards,
Peter

P.S.: the notbook (without comments and packed with 7-zip (see
http://www.7-zip.org/ )) can be found at
http://freenet-homepage.de/Peter_Berlin/Mathe/Interpogration.7z .