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 .
