Re: Algebraic problem solved by simulation
- To: mathgroup at smc.vnet.net
- Subject: [mg53652] Re: [mg53604] Algebraic problem solved by simulation
- From: DrBob <drbob at bigfoot.com>
- Date: Fri, 21 Jan 2005 06:36:44 -0500 (EST)
- References: <BAY1-F248C1FA3357A059D03A43AC1810@phx.gbl>
- Reply-to: drbob at bigfoot.com
- Sender: owner-wri-mathgroup at wolfram.com
Here are the mean AND standard deviation:
a = Rationalize at {0.017, 0.013, 0.012};
delta = 1/20;
Timing[
prob = NIntegrate[Boole[1 -
delta < a.{x, y, z} < 1 + delta]Times @@ a, Sequence @@ Transpose@{{x, y,
z}, 0a, 1/a} // Evaluate]; meanX = NIntegrate[(x + y + z)Boole[
1 - delta < a.{x, y, z} < 1 + delta] Times @@
a, Sequence @@ Transpose@{{x, y, z}, 0a, 1/a} // Evaluate]/prob;
meanXX = NIntegrate[(x + y + z)^2Boole[1 - delta < a.{x, y,
z} < 1 + delta] Times @@ a, Sequence @@ Transpose@{{x, y,
z}, 0a, 1/a} // Evaluate]/prob;
{meanX, Sqrt[meanXX - meanX^2]}]
{0.047 Second,{73.145,5.60071}}
Bobby
On Thu, 20 Jan 2005 21:38:31 +0100, Guillermo Sanchez <guillermosanchezleon at hotmail.com> wrote:
> Dear Bod
>
> If you compare both method the Montecarlo is clearly faster. Also, becouse the
> equation has many solutions (in fact infinite) I need also the standard
> deviation
>
> Thanks for your answer
>
> Guillermo Sanchez
>
> Email (home): GuillermoSanchezLeon at hotmail.com
>
> web: http://web.usal.es/~guillerm
>> From: DrBob <drbob at bigfoot.com> >Reply-To: drbob at bigfoot.com >To: drbob at bigfoot.com,
To: mathgroup at smc.vnet.net
> "Guillermo Sanchez" <guillerm at aida.usal.es>,mathgroup at smc.vnet.net >Subject: [mg53652] Re:
> [mg53604] Algebraic problem solved by simulation >Date: Thu, 20 Jan 2005 14:13:30
> -0600 > >Ah! But NIntegrate is faster: > >a = Rationalize at {0.017, 0.013, 0.012};
>> delta = 1/20; >Timing[num = NIntegrate[( >x + y + z)Boole[1 - delta < a.{x, y,
> z} < 1 + delta] > Times @@ a, Sequence @@ Transpose@{{x, y, z}, 0a, 1/a} // >Evaluate];
>> den = Integrate[Boole[1 - delta < a.{x, > y, z} < 1 + delta]Times @@ a,
> Sequence @@ > Transpose@{{x, y, z}, 0a, 1/a} // Evaluate]; > num/den] >{0.578
> Second, 73.145} > >Bobby > >On Thu, 20 Jan 2005 13:58:25 -0600, DrBob <drbob at bigfoot.com>
> wrote: > >>>>{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 > h
--
DrBob at bigfoot.com
www.eclecticdreams.net