Re: Algorithm

*To*: mathgroup at smc.vnet.net*Subject*: [mg30087] Re: [mg30072] Algorithm*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Fri, 27 Jul 2001 03:52:22 -0400 (EDT)*References*: <200107260520.BAA01183@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Yannis.Paraskevopoulos at ubsw.com wrote: > > Hi all, > > lately I've been trying to translate to mathematica the following > algorithm and it seems that everything is horribly slow and frustrating > to a degree that I think that your kind help will be necessary. > > The algorithm goes like that: > > 1)define t_i=2*i/100 ,i=1,2,...,500 > and X_t0=X > > 2)for j=1,2,...,10,000 > > Do > a) generate x_i from N(0.05,1),p_i from N(0,1) and > > let y_i be 0, with probability 0.9 and 1 with probability > 0.1 > b) calculate ln(X_t(i))=ln(X_t(i-1))+(x_i)+(y_i)*(p_i) > > c) for the smallest i<=500 such that ln(X_t(i))<=0 set > W_j=1-0.4*X_t(i). Otherwise W_j=0 > > 3) calculate final=1-Sum[W_j/10,000, {j,1,10,000}] > > any help will be much appreciated. > > best regards > > yannis > [...] I think you may be having problems for a few reasons. First, you may not be familiar with some relevant Mathematica standard add-on package functions that can handle the generation according to distributions. Second, your specification of the problem could stand some refinement. For example, so far as I can tell, the "t" array (time steps?) is really not needed; we can work with a common index. Moreover there is an initial value, "X", not provided. Also the phrase "Otherwise W_j=0" is unclear; I am guessing you meant "...if no such i exists"? I gave an initial value of 0 for the logarithm of what you refer to as X_t0 (because it is primarily the logs we work with anyway). Making whatever assumptions seem reasonable to interpret the algorithm description, I get the code below. It ran in around 12 minutes on a 300 MHz machine, in version 4.1 of Mathematica. In[1]:= <<Statistics` In[2]:= len = 500; In[3]:= biglen = 10000; In[4]:= initlogbigxx = 0; In[5]:= Timing[ ww = Table[ xx = RandomArray[NormalDistribution[.05,1], {len}]; pp = RandomArray[NormalDistribution[0,1], {len}]; yy = RandomArray[BernoulliDistribution[.1], {len}]; prev = initlogbigxx; logbigxx = Table[ prev = prev + xx[[i]] + yy[[i]]*pp[[i]], {i,len}]; smallindx = Position[logbigxx, _?Negative, 1, 1]; If [smallindx==={}, 0, 1 - 0.4*Exp[logbigxx[[smallindx[[1,1]]]]]], {biglen}]; ] Out[5]= {716.35 Second, Null} In[7]:= InputForm[Take[ww,20]] Out[7]//InputForm= {0.9179154777814722, 0.8966881389164947, 0.9604955637123224, 0.6978081187113281, 0.6387668067441197, 0.9383270710223752, 0, 0.8077776474444456, 0.9758281835484863, 0.6906528869276962, 0.8691539653461342, 0.7665120788037165, 0.9481427893555122, 0, 0.6751854977258047, 0.6533416011239968, 0.70666992443381, 0.9249455563234696, 0.6909018312677464, 0.7051942355651819} This 12 minute run computed 2000 values each time through a loop of 10000 iterations. There may be ways to make it faster (e.g. if the first i with log(Xi) tends to be small, generate xx et al values one at a time until we find our negative value). That said, considering that alot of the computations involve generation of random values, this is fairly good speed. One sees the occasional response to the effect that Mathematica cannot do effective numerical simulations and dedicated numerical software must instead be used. For a problem such as the one above I would not take such a suggestion very seriously. Daniel Lichtblau Wolfram Research

**References**:**Algorithm***From:*Yannis.Paraskevopoulos@ubsw.com