4
$\begingroup$

I'm trying to get a random solution for an equation with arbitrary number of parameters which has infinite solutions. Note that the answer must be in n dimensions where n is the number of parameters of the function.


My tries:

  • The function can be something like x^2 + y^2 + z^2 == 1, in this case, the solution is a sphere and every point that satisfies that equality is a valid solution. This creates a region that consists of valid solutions, so if we generate a random point in that region we have a valid answer for this problem. That can be done using ybeltukov's answer's RegionDistribution:

    implregion = ImplicitRegion[x^2 + y^2 + z^2 + t^2 == 1, {x, y, z, t}];
    region = DiscretizeRegion[implregion];
    pts = RandomVariate[RegionDistribution[region], 100]
    

    but this raises an error saying DiscretizeRegion only takes arguments with less than 4 dimensions.

  • I tried using NSolve to get all the points:

    NSolve[x^2 + y^2 + z^2 == 1, {x,y,z}] (* I don't know what I was supposed to see here *)
    

    but this gives a warning message that says:

     Warning: NSolve::infsolns: "Infinite solution set has dimension at least 1.
      Returning intersection of solutions with (171802*x-113492*y-121484*z)/178835 == 1."
    

    then it outputs some complex valued solutions, which are not valid (in this context):

    {{x->0.252383 - 0.412144I, y->-1.17075 - 0.0795353I, z->-0.0214382 - 0.508549I},
     {x->0.252383 + 0.412144I, y->-1.17075 + 0.0795353I, z->-0.0214382 + 0.508549I}}
    

    As you can see, it gets those wanted infinite solutions, and it even selects a random solution, but it is always the same, and contains complex numbers when what I want is one random answer of that solution set which consists only of real values.

  • I also tried with Reduce doing:

    Reduce[x^2 + y^2 + z^2 == 1]
    

    The result is:

    (x == -Sqrt[1 - y^2 - z^2] || x == Sqrt[1 - y^2 - z^2])
    

    and generating some points doing something like

    Table[
        {RandomChoice[{-1,1}] Sqrt[1 - y^2 - z^2], y, z}/.
          {y-> RandomReal[{-1,1}],
           z-> RandomReal[{-1,1}]},
    {1000}]
    

    gives me incorrect answers because x, y and z are not calculated "at the same time".

Edit:

  • Using FindInstance distributes non uniformly the points, I used

    FindInstance[x^2 + y^2 + z^2 == 1, {x, y, z}, Reals, 1000]
    

    and the result was

    enter image description here

  • And the last method distributes non uniformly too, doing

    Select[Table[{RandomChoice[{-1, 1}] Sqrt[1 - y^2 - z^2], y, z} /.
                     {y -> RandomReal[{-1, 1}], 
                      z -> RandomReal[{-1, 1}]},
           {1000}],
    Im@#[[1]] == 0 &]
    

    gives me

    enter image description here

    this is because the values are calculated depending on the $yz$ plane and in some places x changes faster than in another places.

$\endgroup$
  • $\begingroup$ Have you considered FindInstance? $\endgroup$ – mikado Apr 15 '17 at 18:58
  • 1
    $\begingroup$ The results of your last approach, using Table, certainly are correct. What is your concern with them? $\endgroup$ – bbgodfrey Apr 15 '17 at 19:00
  • 2
    $\begingroup$ Would RandomPoint do what you're after? $\endgroup$ – Chip Hurst Apr 15 '17 at 19:53
  • $\begingroup$ @ChipHurst Yes! $\endgroup$ – Garmekain Apr 15 '17 at 21:29
  • 2
    $\begingroup$ Generating a uniform distribution even on the sphere is highly nontrivial, an ellipsoid is very hard, and you can forget about arbitrary hypersurfaces, so your question is much too general. $\endgroup$ – Igor Rivin Apr 15 '17 at 23:43
2
$\begingroup$

FindInstance works quickly but seems to biased toward points that lie in the xy-plane.

SeedRandom[42]; 
With[{n = 10}, 
  (FindInstance[x^2 + y^2 + z^2 == 1, {x, y, z}, Reals, n] // N)[[All, All, 2]]]
{{0.790419, -0.108434, 0.602893}, {0.367265, 0.930116, 0.}, 
 {-0.538922, -0.842356, 0.}, {-0.538922, 0.842356, 0.}, 
 {0.55489, -0.831924, 0.}, {0.968064, 0.158416, -0.194311}, 
 {-0.407186, -0.913345, 0.}, {0.978044, 0.0658436, 0.197724}, 
 {0.946108, 0.323852, 0.}, {0.55489, 0.831924, 0.}}

So I recommend

SeedRandom[42];
With[{n = 10},
  Module[{y, z},
   z = RandomReal[{-1, 1}, n];
   y = RandomReal[{-1 + #, 1 - #}] & /@ Abs[z];
   MapThread[{Sqrt[1 - #1^2 - #2^2], #1, #2} &, {y, z}]]]
{{0.899547, -0.41092, -0.148189}, {0.972747, 0.0791143, -0.217954}, 
 {0.903025, 0.301654, -0.305861}, {0.882229, 0.46164, -0.0925187}, 
 {0.760152, 0.640032, 0.111927}, {0.701265, 0.574829, -0.421661}, 
 {0.868464, 0.284055, -0.406304}, {0.714105, -0.38114, -0.587185}, 
 {0.922918, 0.161119, -0.349661}, {0.321531, -0.0217459, 0.946649}}
$\endgroup$
  • $\begingroup$ It is also biased, two peaks seem to have more points than the other two; check it out $\endgroup$ – Garmekain Apr 15 '17 at 21:22
1
$\begingroup$

Region functionality in M11.2+ can handle Ellipsoid objects. For example:

pts = RandomPoint[
    RegionBoundary @ Ellipsoid[{0,0,0,0},{1,2,3,4}],
    10
]

{{0.702264, -0.506915, -1.00603, -2.29828}, {0.672747, 1.1246, 0.948138, -1.44966}, {-0.543024, 0.140708, -2.00356, -2.01652}, {0.652806, -1.37052, -0.93464, -0.339369}, {-0.38219, -0.873541, -2.1748, 1.48397}, {0.220089, -1.55027, 1.56574, -1.11951}, {-0.38265, 0.746885, 2.49541, -0.596283}, {0.52443, 1.36327, -1.48924, 0.471912}, {-0.212721, 1.43821, 1.50939, 1.71813}, {0.649815, -0.909925, 0.557857, -2.31921}}

Check that they all lie on the Ellipsoid:

Norm/@(pts.DiagonalMatrix[1/Range[4]])

{1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.