7
$\begingroup$
 a = RandomVariate[UniformDistribution[{0, 1}]];
 b = RandomVariate[UniformDistribution[{0, 1}]];
 c = RandomVariate[UniformDistribution[{0, 1}]];
 a + a b + a b c

I want to continue picking random reals in (0,1) and adding the product to the previous sum. How can Mathematica help me here?

Best regards Geoffrey Critzer

$\endgroup$
13
$\begingroup$

Here's one way:

SeedRandom[1];
a = RandomVariate[UniformDistribution[{0, 1}]];
b = RandomVariate[UniformDistribution[{0, 1}]];
c = RandomVariate[UniformDistribution[{0, 1}]];
a + a b + a b c

0.980367

SeedRandom[1];
values = RandomVariate[UniformDistribution[{0, 1}], 3];
Total@FoldList[Times, values]

0.980367

The number 3 can be replaced by any number, however many times you want to iterate.

Here's a procedural solution (with the definition of values as in the previous example):

prod = First[values];
sum = First[values];
Do[
  prod *= v;
  sum += prod,
  {v, Rest[values]}
  ];
sum

0.980367

$\endgroup$
8
$\begingroup$

C.E.'s answer is great already. I would just like to point out that we may exploit here that floating point addition is usually significantly faster than floating point multiplication that FoldList is just slow, and that multiplication can be cast into addition by applying Log so that we can use Accumulate instead. Morever, we may use vectorized built-in routines for that.

n = 1000000;
values = RandomVariate[UniformDistribution[{0, 1}], n];

r1 = Total@FoldList[Times, values]; // RepeatedTiming // First
r2 = Total[Exp[Clip[Accumulate[Log[values]], {-700., ∞}]]]; // RepeatedTiming // First

Max[Abs[r1 - r2]]

0.070

0.0053

0.

For those who wonder what the Clip is for: This is in order to prevent underflow error handling to occurr (the latter slows down things considerably); that happens at about Exp[-709.] or so.

Edit

It is even faster to write a short compiled version of C.E.'s procedure (if do not count in the compilation time):

cf = Compile[{{x, _Real, 1}},
   Block[{prod = 1., sum = 0.},
    Do[prod *= Compile`GetElement[x, i]; sum += r, {i, 1, Length[x]}];
    sum
    ],
   CompilationTarget -> "C"
   ];

Now:

r3 = cf[values]; // RepeatedTiming // First
Max[Abs[r1 - r3]]

0.0013

1.77636*10^-15

Remark

I formerly claimed that floating point multiplication were slower than floating point addition. As Roman pointed out, that is not correct. While multiplication probably has higher complexity (and with floating point computations, some quite counterintuitive things happen), modern hardware is built such that variuous steps of the multiplication are performed in parallel. Nowadays, there is even a single circuit for fused multiply-add (FMA) and not necessarily any separated addition circuit, so addition and multiplication should take basically the same time.

$\endgroup$
  • 2
    $\begingroup$ On a modern CPU, multiplication is not slower than addition. See, e,g., this discussion that's already ten years old. Single-clock-cycle multipliers have been around for a while now. $\endgroup$ – Roman Aug 10 at 19:04
  • $\begingroup$ @Roman Yes, actually, I've read many things like that today, in particular about things that were introduce with the Haswell. Now I am not sure anymore why the Accumulate code is faster than the FoldList. $\endgroup$ – Henrik Schumacher Aug 10 at 19:40
  • $\begingroup$ @Roman I am still working with a Haswell processor and there the situation is such that both ADD and FMA (which is used there for multiplication on this architecture IIRC; since Skylake, FMA is used for both) have the same throughput but ADD has significantly less latency. And that might make the difference: We are not multipling or adding two vectors componentwise (in which case it is only about throughput, not about latency), we are indeed folding here - and that should also depend on latency, IMHO. $\endgroup$ – Henrik Schumacher Aug 10 at 19:43
  • $\begingroup$ Maybe the sole reason for the timing differences here is the simple fact that there is no such things as an efficiently implemented multiplicative variant of Accumulate in Mathematica... $\endgroup$ – Henrik Schumacher Aug 10 at 19:44
  • $\begingroup$ Maybe related to this discussion about the memory hogging of NestWhile? Could it be that in general, nesting and folding are just not that efficient? After all, Accumulate[x] is ten or twenty times faster than FoldList[Plus, x], whereas FoldList[Times, x] is only imperceptibly slower than FoldList[Plus, x]. $\endgroup$ – Roman Aug 10 at 19:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.