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
- References:
- Algebraic problem solved by simulation
- From: guillerm@aida.usal.es (Guillermo Sanchez)
- Algebraic problem solved by simulation