       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:=
<< "Imtek`Interpolation`"
<< "NumericalMath`NIntegrateInterpolatingFunct`"

Here we get a disturbed array of points:

In:= SeedRandom;
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:=
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:=
{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=
{0.33819365081907227, 0.5860641465257586}
Out=
1.1739273460216542

just to visualize the current situation:

In:=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,
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}]}], ImageSize -> 500]]

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

And here we go:

In:=
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=
3.079901416226034

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

percent is quite good, isn't it?

(just to verify, that Pi should be the outcome:
In:=
Block[{r0 = Sqrt[2*Log], 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=
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 .

```

• Prev by Date: Re: asymptotics
• Next by Date: Re: Quadratic form: symbolic transformation
• Previous by thread: ListIntegrate3D ?
• Next by thread: Export issue