Bounds on random number

Question

I have a Markov Chain Monte Carlo code that uses a chi squared function to fit some data. I am running into an issue now with getting unrealistic values for some of the free parameters in the code. What I would like to do is place an lower and upper bound on the parameter value (say k1[[i]] below ). Each iteration (k1[[i]]) of the function produces a random number that is added or subtracted from k1[[i-1]] giving the next step in the random walk. I would like to be able to test each current parameter (k1[[i]]) and if it is above or below a cutoff value, choose another random number until I return to the specified range.

If k1[[i]] < lowest or k1[[i]] > largest, choose another Random number until lowest < k1 < largest

However I am not quite sure how to implement this into the code for each parameter.
Does anyone have any suggestions?

Table[
  {

   k1[[i]] =  k1[[i - 1]] + weight*k1[[i - 1]]*RandomReal[{-1.0, 1.0}],

   k2[[i]] = k2[[i - 1]] + weight*k2[[i - 1]]*RandomReal[{-1.0, 1.0}],

   k3[[i]] = k3[[i - 1]] + weight*k3[[i - 1]]*RandomReal[{-1.0, 1.0}]
   },

 {i,1,range}

]

Can you use a While loop and get it to choose random numbers until the condition in the while loop (the number being out of the desired bound) is no longer met. — Jonathan Shock, May 16 '13 at 00:00
How about setting the limits of RandomReal such that they will only return an allowable value? If you know the values of lowest and largest then you should be able to get your limits from the inequalities. — bobthechemist, May 16 '13 at 00:03
@bobthechemist, I'm a bit confused because user7163 seems to know about these bounds as they are being used in the code example... — Jonathan Shock, May 16 '13 at 00:11
Hi, I believe the while loop should work. I'll try that. @bobthechemist I do not know the range of the parameter values, however I was getting values smaller than the smallest number available on mathematica. — tarhawk, May 16 '13 at 02:35

Mr.Wizard · Answer 1 · 2013-07-30T20:11:21.857

Your three lists k1, k2, k3 are independent and constructed similarly so I shall show how to create one of them.

Since you will be reusing a result from each prior calculation it is natural to use NestList.
I will write a makeNext function which takes two arguments:

a weight value
a test function which the output must pass

And generates a new function which is to be given the value of the last valid result (i.e. k[i-1]).

makeNext[weight_, test_] :=
 Compile[{k0},
  Module[{x = 0`},
   While[ ! test[ x = k0 + weight k0 RandomReal[{-1, 1}] ] ];
   x
  ]
 ]

Applying this function with a weight of 0.7 and a bound of (0.3, 0.4) to a starting value of 1`:

NestList[makeNext[0.7, 0.3 < # < 0.4 &], 1`, 15]

{1., 0.331259, 0.369334, 0.34816, 0.337829, 0.360368, 0.34552, 0.379429, 0.330914,
 0.396059, 0.359248, 0.311805, 0.335917, 0.313851, 0.319714, 0.32371}

You can use Rest to drop the first value if it is unwanted.

score 3 · Accepted Answer · answered Jul 30 '13 at 21:36

Note For the case of a uniform random distribution, it will be equivalent to precompute the bounds on the random variable:

makeNext[weight_,xl_,xh_] := Compile[{k0}, Module[{}, 
       k0 + weight k0 RandomReal[{
           Max[-1, (xl - k0)/(weight k0 )],
           Min[1, (xh - k0)/(weight k0 )]}]]]
NestList[makeNext[0.7,.3,.4], 1`, 15]

If you have a nonuniform distribution you'd have a bit of work to constrain the range while maintaining the same statistics.

Just for fun maybe you want to bounce off the boundary..

makeNext[weight_, xl_, xh_] := Compile[{k0}, Module[{x}, 
      x = k0 + weight k0 RandomReal[{-1, 1}];
      While[ x > xh || x < xl,
      x = Which[x > xh, 2 xh - x, x < xl, 2 xl - x]]; x]]

Good observation. I assumed the example was not entirely representative and I tried to keep my answer general, but this is a lot better for the specific case. It reminds me of method 3 in my own answer to a different problem. — Mr.Wizard, Jul 30 '13 at 21:47

Bounds on random number

2 Answers2