# How to improve this function that selects parameters from arbitrary distribution using Recursion?

Below, I have defined a function using Which to select model parameters based upon various conditions. My questions are: Is there any danger in writing this function using Recursion? Is there a more efficient way to write this function that does not require recursion? I'm also only showing an example Uniform distribution here, however the distributions used in the actual program will be a derived distribution from some input data.
1. The first condition selects random variates from a distribution.
2. The second condition selects random variates from a different distribution.
3. The third condition performs a random walk on the random variate selected from 2.

The function below is called repeatedly and after each call I determine if the algorithm has selected a better set of parameters given some fitting criteria. This particular model fitting strategy is known as Sequential Monte Carlo or Particle Filtering. While I have been using standard MCMC, my goal is now to maximize the coverage of parameter space in an efficient manner.

Clear[selectParameters,startingDistributions,y]
m2={k1p,km1};
startingDistributions[model_]:={UniformDistribution[{0.,50.}],UniformDistribution[{0.5,100.}]};
selectParameters[model_,epsilon_Integer,i_]:=
y[i]=
Which[
epsilon==1,
(model[[#1]]->RandomVariate[startingDistributions[model][[#1]]]&)/@Range[Length[model]],
epsilon>1&&i==1,
(model[[#1]]->RandomVariate[startingDistributions[model][[#1]]]&)/@Range[Length[model]],
epsilon>1&&i>1,
(model[[#1]]->y[i-1][[#1,2]]+0.025 y[i-1][[#1,2]] RandomReal[{-1,1}]&)/@Range[Length[model]]]
test=Table[selectParameters[m2,2,i],{i,1,100}];
lp={k1p,km1}/.test;
ListPlot[{lp[[All,1]],lp[[All,2]]},Joined->True]


Nothing intrinsically wrong with what you've done, but perhaps a bit messy. For example, here's the epsilon>1 case done a bit more Mathematica style. Easier to read and faster:

(* epsilon>1 case as function *)

genResults = Module[{inc = 1, init = RandomVariate /@ #1,
rvs = RandomReal[{-1, 1}, {#2, Length@#1}]},
NestList[# + .025 # rvs[[inc++]] &, init, #2 - 1]] &;

(* The above in use: make a dist list, and how many iterations... *)

sd = {UniformDistribution[{0., 50.}], UniformDistribution[{0.5, 100.}], UniformDistribution[{10., 75.}]};

n = 10000;

(* Generate the walk and plot *)
result = genResults[sd, n];

ListPlot[Transpose@result, Joined -> True]


If you run yours after doing say SeedRandom[1], then the above after seeding the same, you'll see that your lp and the results from above are exactly the same.

Per your comments, here's an alternative way to replace the whole epsilon>1 case with an easier to read code chunk that is perhaps easier for you to add arbitrary actions:

sd = {UniformDistribution[{0., 50.}], UniformDistribution[{0.5, 100.}]}

f[1] := RandomVariate /@ sd;
f[n_] := f[n] = f[n - 1] + .025 f[n - 1] RandomReal[{-1, 1}, Length@sd];


You can then generate the results for some number of steps (100 here) via:

f/@Range@*100


and the various f will remain defined to do with as you wish (just don't forget to clear them when you're done with them!)

As I said, nothing wrong with your way, use whatever style floats your boat!

• Hi @rasher, Thanks for the response. If you have time, could you explain the pure function approach in a little more detail. This is something I am trying to master as well. Also, suppose I wanted to compute something after each selection of random variates, how could I do that with your method above? For example suppose as a test, I wanted to too just Sum the iterates and test them against the previous iterate? Total[[i]] vs. Total[[i-1]]. The method I have written above allows this, but as you mention it is messy and inefficient. Feb 22, 2015 at 13:34
• @tarhawk: I'm not sure I get what you're asking - the prior iterate is the slot (#) in the function, the new one is the result of the function contents, so you could just jigger the function to save the prior, generate the new, then do whatever comparisons/tests you need. Perhaps update OP with detail on what you want to accomplish? That said, I'll update the post with an alternative and still cleaner way that retains the indexed variable construct you use.
– ciao
Feb 22, 2015 at 21:45
• I have updated the OP and provided an example of the output from a few calls to my original function. Perhaps this will help clarify my goal. I note that this function works, however as you've shown above, I would like to learn how to more efficient write such functions using a pure function style. Feb 23, 2015 at 14:40