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
- Algorithm