MathGroup Archive 2005

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

Search the Archive

Re: Algebraic problem solved by simulation

  • To: mathgroup at smc.vnet.net
  • Subject: [mg53646] Re: [mg53604] Algebraic problem solved by simulation
  • From: DrBob <drbob at bigfoot.com>
  • Date: Fri, 21 Jan 2005 06:36:12 -0500 (EST)
  • References: <200501200847.DAA04089@smc.vnet.net>
  • Reply-to: drbob at bigfoot.com
  • Sender: owner-wri-mathgroup at wolfram.com

>> {X1, ...Xn} follows a Uniform random  distribution

You've drawn X1, X2, and X3 from THREE uniform distributions, with the same lower limit but three different upper bounds. Was that your intent? If so, and if you mean to compute the CONDITIONAL mean of x+y+z given that a.{x,y,z} is between the two bounds, then...

Here's an exact method that closely agrees with random results from your code:

a = Rationalize[{0.017, 0.013, 0.012}];
delta = 1/20;
num = Integrate[(x + y + z)*Boole[1 - delta <=
        a . {x, y, z} <= 1 + delta]*Times @@ a,
     Sequence @@ Transpose[{{x, y, z}, 0*a, 1/a}]];
den = Integrate[Boole[1 - delta <= a . {x, y, z} <=
        1 + delta]*Times @@ a, Sequence @@
      Transpose[{{x, y, z}, 0*a, 1/a}]];
num/den

23864575/326264

N[%]

73.145

(Boole is in the System context now, but used to be in the Calculus`Integration` package. Load that if you have an older version.)

In both integrals, Times@@a dx dy dz is the volume element of integration. Boole is 1 when the dot product is right, 0 otherwise. In the first integral we compute the integral of x+y+z WHEN the dot product is right, and in the second we compute the probability of that being so.

This may be no faster than your Monte Carlo method, however.

Bobby

On Thu, 20 Jan 2005 03:47:44 -0500 (EST), Guillermo Sanchez <guillerm at aida.usal.es> wrote:

> (* Dear group;
> Given the equation L1 < a1 X1, ..., an Xn < L2 where {a1, ..an} are
> known and {X1, ...Xn} follows a Uniform random  distribution I wish to
> find the average and the estandar deviation of X, being X = X1 + ... +
> Xn.
>
> For instance, for solving 0.95 < 0.017 X1 + 0.013 X2 + 0.012 X3 < 1.05
> I apply this method*)
>
> a = {0.017, 0.013, 0.012};
>
> rr := {Random[Real, {0, 1/0.017}], Random[Real, {0, 1/0.013}],
> Random[Real, {0, 1/0.012}]};
>
> eq1 := (X = rr; If[0.95 < a . X < 1.05, Plus @@ X] ) ;
>
> L = Array[eq1 &, 1000] ;
>
> sol1 = Cases[L, x_?NumericQ];
>
> {Mean[sol1], StandardDeviation[sol1]}
>
> (*Can any body find a better method?, In the real problems the size of
> {X1, ..., Xn} is too big (1000)
>
> Guillermo*)
>
>
>
>



-- 
DrBob at bigfoot.com
www.eclecticdreams.net


  • Prev by Date: Re: ShowLegend....
  • Next by Date: Re: Re: random matrix from row and column sums
  • Previous by thread: Algebraic problem solved by simulation
  • Next by thread: Re: Algebraic problem solved by simulation