# Uniform distribution on z = 1 − √( x ^2 + y ^2) , z ≥ 0 [duplicate]

Let $S$ belong to ${\rm I\!R^3}$, $z = 1 − \sqrt{x^2 + y^2}$ , $z ≥ 0$.

I need to write a program that generates a random point on $S$ with uniform distribution. I've tried to use

g[z_, fi_] = {(1 - z)*Cos[2 Pi*fi], (1 - z)*Sin[2 Pi*fi], z};
sample[n_] := Table[g[RandomReal[], RandomReal[]], {n}];

sampleX[n_] := Table[g[RandomReal[], RandomReal[]][[1]], {n}];
sampleY[n_] := Table[g[RandomReal[], RandomReal[]][[2]], {n}];
sampleZ[n_] := Table[g[RandomReal[], RandomReal[]][[3]], {n}];


But the result isn't exact, because I have more points on the top and less points near $OXY$.

How could I make a Uniform distribution on this surface?

• One more question :) How to calculate an integral(z dS) using wolfram mathematica? Commented Sep 27, 2017 at 9:29
• Welcome! To make the most of Mma.SE start by taking the tour now. It will help us to help you if you write an excellent question. Edit if improvable, show due diligence, give brief context, include minimal working example of code and data in formatted form. As you receive give back, vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. Commented Sep 27, 2017 at 9:42
• Your "one more question" can be found in the documentation of Integrate. Follow the link. Commented Sep 27, 2017 at 9:44
• w.r.t. t = 1-z, so the cumulative distribution function is t^2. Why is it so? Commented Sep 27, 2017 at 12:56

If you set up your region S as a region in Mathematica you can use RandomPoint to generate uniformly distributed points over it.

reg = ImplicitRegion[
z == 1 - Sqrt[(x^2 + y^2)], {{x, -1, 1}, {y, -1, 1}, {z, 0, 1}}];
dreg = DiscretizeRegion[reg]
rp = RandomPoint[dreg, 10000];
Show[ContourPlot3D[
z == 1 - Sqrt[(x^2 + y^2)], {x, -1, 1}, {y, -1, 1}, {z, 0, 1}],
Graphics3D[{PointSize[0.05], Point[rp]}]]


This should be fine for most purposes, but note that the discretised region dreg is a numerical approximation to the ImplicitRegion, reg. So, for example

point = RandomPoint[dreg]
RegionMember[dreg, point]
RegionMember[reg, point]

(* {-0.538245, -0.716972, 0.103292}
True
False *)


That is, while point in the discretised region, dreg, it is not actually in the implicit region reg. Although it is close:

RegionDistance[reg, point]

(* 0.000129349 *)


The approximation of the mesh in dreg could be improved with MaxCellMeasure.

In case you're wondering, sadly we can't just generate RandomPoints from reg:

RandomPoint[reg]

(* RandomPoint[ImplicitRegion[z == 1 - Sqrt[x^2 + y^2] && -1 <= x <= 1 && -1 <= y <= 1 && 0 <= z <= 1, {x, y, z}]] *)


It "should" work (according to documentation), but Region functionality in Mathematica often has funny quirks to it.

UPDATE: Comparison of approaches to uniform distributions over a region embedded in a higher dimensional space

There have been two general approaches to generating uniformly distributed random points over the surface of a cone. @Coolwater and @george2079 opted for analytical solutions, while @Carl Woll and myself used Mathematica's built in Region functionality. I wanted to have a look at how these approaches compare in accuracy and speed; and get a sense of how best to approach similar problems.

An important feature of this question is that the region dimension is less than the embedding dimension (and this is the reason for the woeful accuracy of RandomPoint applied to DiscretizeRegion; if RegionDimension == EmbeddingDimension accuracy shouldn't be a problem). With this in mind, the take-home message based on the table below seems to be:

1. Keep to analytical solutions as far as possible. A little math goes a long way.
2. If you must use Regions, use the ParametricRegion or built-in ones.
3. Use ImplicitRegion and DiscretizeRegion only in emergencies.

First, a few definitions:

carlreg =
RegionIntersection[RegionBoundary[Cone[{{0, 0, -1}, {0, 0, 1}}, 2]],
ConicHullRegion[{0, 0, 0}, {{1, 0, 0}, {0, 1, 0}}, {{0, 0, 1}}]];
implicitreg =
ImplicitRegion[
z == 1 - Sqrt[x^2 + y^2], {{x, -1, 1}, {y, -1, 1}, {z, 0, 1}}];
discreg = DiscretizeRegion[implicitreg];
parareg =
ParametricRegion[{Cos[θ] u, Sin[θ] u, 1 - u},
{{θ, 0, 2 π}, {u, 0, 1}}];

coolwater[n_] :=
Transpose[{Cos[#] #2, Sin[#] #2, 1 - #2} &[RandomReal[2 π, n],
Sqrt[RandomReal[1, n]]]]
george[n_] :=
Append[#1 Through[{Sin, Cos}@#2], 1 - #1] & @@@
Transpose[{RandomVariate[ProbabilityDistribution[2 r, {r, 0, 1}],
n], RandomReal[{0, 2 π}, n]}]
carlwoll[n_] := RandomPoint[carlreg, n]
aardvark[n_] := RandomPoint[discreg, n]
parametric[n_] := RandomPoint[parareg, n]


Now generate some points

{names, timing, points} =
Transpose[Join[{#}, AbsoluteTiming[#[1000]]] & /@
{coolwater, george, carlwoll, aardvark, parametric}];
regdist = Mean[RegionDistance[implicitreg, #]] & /@ points;
solerror = Mean[Abs[#3 - (1 - Sqrt[#1^2 + #2^2])] & @@@ #] & /@ points;
TextGrid[Join[{{"", "Timing", "RegionDistance", "Solution error"}},
Transpose[{names, timing, regdist, solerror}]]]


\begin{array}{|c|ccc|} \hline \text{} & \text{Timing} & \text{RegionDistance} & \text{Solution error} \\ \hline \text{coolwater} & 0.0003290 & 1.25468 \times 10^{-10} & 1.68754 \times 10^{-17} \\ \text{george} & 0.0049363 & 1.66830 \times 10^{-10} & 1.87628 \times 10^{-17} \\ \hline \text{carlwoll} & 0.0458737 & 5.67581 \times 10^{-10} & 6.13953 \times 10^{-17} \\ \text{aardvark} & 0.0006388 & 0.000600271 & 0.000849862 \\ \text{parametric} & 0.0251807 & 7.21428 \times 10^{-11} & 1.74305 \times 10^{-17} \\ \hline \end{array}

A few words about the difference between implicitreg, parareg and @CarlWoll's carlreg. They are clearly the same mathematical region, their Region-based properties (RegionMeasure, RegionCentroid, etc) are identical, and FindInstance cannot find any points that they do not all have in common. However, RandomPoint can be used directly on carlreg and parareg but not on implicitreg. I don't have an explanation why that is the case, but it is possibly related to this question by @rhermans. Basically, ImplicitRegions are much more complicated for Mathematica to work with. So using them for a case like this, where there are many other options, is clearly (in hindsight) a mistake.

• I think that problem is here: g[z_, fi_] = {(1 - z)*Cos[2 Pifi], (1 - z)*Sin[2 Pifi], z}; I didn't use Jacobian to generate the points, but idk how to do that Commented Sep 27, 2017 at 8:44
• Is it real to get the coordinates of these points, which are generated? Commented Sep 27, 2017 at 8:48
• Ahh, ofc! Just there is a list rp Commented Sep 27, 2017 at 8:50
• Thanks a lot! Works Commented Sep 27, 2017 at 8:51

The radius increase linearly w.r.t. t = 1-z, so the cumulative distribution function is t^2. Hence transforming uniform variables for t by the square root (the quantile) makes the result uniform.

samp[n_] := Transpose[{Cos[#] #2, Sin[#] #2, 1 - #2} &[
RandomReal[2 π, n], Sqrt[RandomReal[1, n]]]]

surf = ContourPlot3D[z == 1 - Sqrt[x^2 + y^2], {x, -1, 1}, {y, -1, 1}, {z, 0, 1}];
Show[surf, ListPointPlot3D[samp[9000]]]


• w.r.t. t = 1-z, so the cumulative distribution function is t^2. Why is it so? Commented Sep 27, 2017 at 12:08
• The radius of the circles around the z-axis is t. And the radius is proportional to the density of points conditional on z (or t same-same). The CDF t^2 is the result of integrating the radius (wrongly scaled density) and scaling to have range 0 to 1. Commented Sep 27, 2017 at 13:14
• Thanks! Understood Commented Sep 27, 2017 at 14:05

Your surface is just the boundary of a cone (excluding the base), and RandomPoint will work with the region boundary of a Cone object. We can intersect the region boundary of a slightly extended cone with a region with $z>0$ to eliminate the base. RandomPoint will also work with this region intersection object. So, I will use a cone with a circle of radius 2 at $(0, 0, -1)$ as the base and $(0, 0, 1)$ as the top, so that the intersection with $z>0$ is the desired cone:

reg = RegionIntersection[
RegionBoundary[Cone[{{0, 0, -1}, {0, 0, 1}}, 2]],
ConicHullRegion[{0, 0, 0}, {{1, 0, 0}, {0, 1, 0}}, {{0, 0, 1}}]
];

pts = RandomPoint[reg, 1000];

Graphics3D[{Point@pts, RGBColor[0, 1, 1, .5], Cone[{{0, 0, 0}, {0, 0, 1}}]}]


And, here is a check that pts satisfy the desired equation:

Tally @ Chop[Function[{x, y, z}, z - (1 - Sqrt[x^2+y^2])] @@@ pts]


{{0, 1000}}

small tweak using RandomVariate :

n = 10^3
ListPointPlot3D[
Append[#1 Through[{Sin,Cos}@#2], 1 - #1] & @@@
Transpose[{
RandomVariate[ProbabilityDistribution[2 r , {r, 0, 1}], n],
RandomReal[{0, 2 Pi}, n]}]]