4
$\begingroup$

I would need some (a lot of) help to improve/rewrite some code.

Since my real problem is a bit too complicated to explain, I'll try to show a simplified scheme of what I want to achieve.

I create some random numbers on an interval between 0 and 1 with:

maximumRandomNumberFirstStep=10000; randomNumberFirstStep=RandomReal[{0,1},maximumRandomNumberFirstStep];

Out of some reasons, I am only interested in values which are higher than a special limit value. So, I am searching for these values with:

limitValue=0.5; helpListFirstStep=ParallelTable[If[randomNumberFirstStep[[i]]<limitValue,i,0], {i,1,maximumRandomNumberFirstStep}];

For all values which are lower than this limit value, I want to generate some new random numbers. At the moment I am doing this with:

maximumRandomNumberSecondStep=maximumRandomNumberFirstStep-Count[helpListFirstStep,0]; randomNumberSecondStep=RandomReal[{0,1},maximumRandomNumberSecondStep];

Now, I want to add the new random numbers (randomNumberSecondStep) to all the old random numbers (randomNumberFirstStep) which are below my limit value. I am doing this with this very ugly code:

j=1; randomNumberSecondStep2=Table[If[ helpListFirstStep[[i]]==0, {randomNumberFirstStep[[i]],0}, {randomNumberFirstStep[[i]]+randomNumberSecondStep[[j]],j++} ],{i,1,maximumRandomNumberFirstStep}]; randomNumberSecondStep2=randomNumberSecondStep2[[All,1]];

(by the way, this code does only work with Table but not with ParallelTable if I execute it more than once; perhaps somebody has also a solution for that)

If one creates a histogram Histogram[randomNumberSecondStep2], it can be seen,= that there are still values which are below my limit value. So, I have to copy my code from above, rename Second to Third and execute it again. I am doing this so often until none of my randomNumberNthStep2-values are below my limit value.

This works, but I think it cannot be programmed worse than I am doing it (copying, renaming etc. is really ugly, especially if your real problem is more complicated because of more variables).

So, how can I write a better program in which all the copying and renaming steps can be avoided? I think there should be a solution in which all these steps are executed automatically. Unfortunately, my programing skills are too low to develop a nice solution.

Improve this question
$\endgroup$
8
  • $\begingroup$ Why doesn't this work: RandomReal[{lowerlimit,1},maximumRandomNumberFirstStep] where lowerlimit is the lower bound you wish to use? $\endgroup$ Commented Mar 30, 2012 at 15:31
  • $\begingroup$ Thanks for the reply. I feared that this would be one of the first answers. That is something I cannot do in my real problem. Simply said, I need the randomNumberNthStep2-values and I need their distribution. $\endgroup$ Commented Mar 30, 2012 at 15:39
  • $\begingroup$ I thought of this at first, but if you look closely, there's some funny additions going on inside the Ifs. It's hidden in the default view and you'll have to scroll to the side to see {randomNumberFirstStep[[i]]+randomNumberSecondStep[[j]],j++}. A NestWhile would probably be of help here. Not going to work on it now though. $\endgroup$ Commented Mar 30, 2012 at 15:41
  • $\begingroup$ @partial81 as pointed out by R.M, I apparently did not read it closely enough. Do you have an idea of what distribution your looking for? If you do, we may be able to bypass the multiple runs. $\endgroup$ Commented Mar 30, 2012 at 15:51
  • $\begingroup$ Thanks for editing my question @R.M. Yes, the addition on the If is important since I need the randomNumberNthStep2-values. And thanks for the hint with NestWhile. If I find a solution by myself (that would be a miracle ;-), I am going to post it for sure! $\endgroup$ Commented Mar 30, 2012 at 15:51

2 Answers 2

Reset to default
8
$\begingroup$

This is wasteful but, if I understand the problem correctly, it does the trick with minimal code using ReplaceRepeated.

n = 10000; limitValue = 0.5; RandomReal[{0, 1}, n] //. x_ /; x < limitValue :> x + RandomReal[{0, 1}]

I was actually surprised at the speed of it.

Edit:

Per request in the comments. To keep track of the number of additions we can modify things slightly with Reap and Sow.

n = 10000; limitValue = 0.5; {res, count} = Reap[Transpose[{RandomReal[{0, 1}, n], Range[n]}] //. {x_, i_} /; x < limitValue :> (Sow[i]; {x + RandomReal[{0, 1}], i})];

Now res[[All,1]] contains the values you want and Tally[count[[1]]] contains a list of pairs {i,k} where i is the index and k is the number of additions.

Improve this answer
$\endgroup$
6
  • $\begingroup$ The only thing is, I think his code does it only once. Perhaps changing //. to /.? $\endgroup$ Commented Mar 30, 2012 at 17:22
  • $\begingroup$ @Rojo I was under the impression that the issue is a need to repeatedly apply the procedure in an automated way. $\endgroup$ Commented Mar 30, 2012 at 17:29
  • $\begingroup$ WOW, that seems to do the thing I want to have! This is really a very smart solution; I think I should be able to adapt it to my main program. Then I will let you and @rcollyer know if it works there too. I think it could be useful to know how often values below the limit value were added up (e.g. is the value 1.4 the sum of 2, 3,... add ups). Do you see also a simple solution for that? I unfortunately do not at this moment. $\endgroup$ Commented Mar 30, 2012 at 19:54
  • $\begingroup$ Dear Andy Ross, Thanks a lot for editing your question! That is really helpful! I just want to get sure if I understand and do everything correctly: If I use Tally[count[[1]], then I get something like: {{1, 1}, {3, 1},...}. This means that at index 2 no addition took place, right? If I want to use my random number inside of a function (lets say Sin[x]+Cos[x]), and the result of this function has to be above the limit value, then I could do something like (line 1,3,4 as in your code) limitValue=1.2; Sin[x]+Cos[x]<limitValue...; resFinal=Sin[res[[All,1]]]+Cos[res[[All,1]]]; Do you agree? $\endgroup$ Commented Mar 31, 2012 at 15:09
  • $\begingroup$ Yes, you interpreted Tally[count[[1]]] correctly. As for the function bit I would define a function f[x_]:= Sin[x]+Cos[x] then you would use {res, count} = Reap[Transpose[{f[RandomReal[{0, 1}, n]], Range[n]}] //. {x_, i_} /; x < limitValue :> (Sow[i]; {x + f[RandomReal[{0, 1}]], i})]. You may need to experiment since x + f[RandomReal[{0, 1}]] and f[x + RandomReal[{0,1}]] give completely different distributions and I don't know which you are wanting. $\endgroup$ Commented Mar 31, 2012 at 16:26
7
$\begingroup$

Here's another method that's still wasteful in the sense of generating more random numbers and performing more arithmetic operations than strictly necessary, but nonetheless substantially faster than the ReplaceRepeated approach (by about a factor of 15). Whether this extra performance is of any use to you I don't know (perhaps, given that you were using ParallelTable to start with):

n = 100000; lowerLimit = 0.5; FixedPoint[ # + RandomReal[{0, 1}, n] * UnitStep[lowerLimit - #] &, RandomReal[{0, 1}, n] ]; // Timing

On my computer this produces a timing of about 0.1 seconds versus around 1.5 seconds for Andy Ross's method (without Sow/Reap).

Edit

If we want to keep track of the number of additions performed in each place, a minor modification can be used. Here it's preferable not to use Sow/Reap, to avoid incurring the overhead of storing the tallies at each individual step, when in fact only the total is needed. As a result the performance of this version is only about 5% worse than that without keeping a tally, which is probably an insignificant difference.

Module[{tally, count = ConstantArray[0, n]}, tally = (count += #; #) &; FixedPoint[ # + RandomReal[{0, 1}, n] * tally@UnitStep[lowerLimit - #] &, RandomReal[{0, 1}, n] ] ~List~ count ]
Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for posting this nice solution although I have already accepted Andy’s answer. You are right, I need a bit extra performance in my program to avoid long calculation times (as you recognized by the use of e.g. ParallelTable). So I could imagine using also your solution. I have to see how good I can implement it in my main program. Thanks again! $\endgroup$ Commented Apr 3, 2012 at 6:29

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.