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