(It might be useful to refer to page 9 of Multi-Armed Bandit Allocation Indices by Gittins, Glazebrook and Weber if you have it, because there explanation will be much better than mine.)I'm trying to write a program that finds that optimal expected reward for a 2-armed Bernoulli bandit process where the first arm has prior probability theta ~ Beta(alpha,beta) of success and the second has a fixed probability p of success, and we have a finite horizon n. The prior is updated via Bayes' theorem after each pull. The reward at the kth pull is gamma^k for a success, and 0 for a failure, where gamma is in [0,1). Using the theory from the book, I've written the following program in R, which utilises dynamic programming:reward = function(alpha,beta,p,gamma,n){if (n==1)    {    return(max(alpha/(alpha+beta),p))    }else    {    success = alpha/(alpha+beta) * (1 + gamma*reward(alpha+1,beta,p,gamma,n-1))    failure = beta/(alpha+beta) * (gamma*reward(alpha,beta+1,p,gamma,n-1))    other = p*(1-gamma^n)/(1-gamma)    return(max(success+failure,other))    }}This program tells you the expected optimal reward, and I'm pretty sure it works, since for small N it agrees with my hand calculations. However, this program does not output the sequence of pulls required to obtain the optimal reward. How should I go about modifying this program so as to produce the list of pulls as well?(For what it's worth, I only know a little of R and MatLab. And sorry if this would be more suited to math.stackexchange)Edit:I had already tried to do as Wandering Logic said in matlab, using the following:function [answer] = reward(alpha,beta,p,gamma,n)global arm karm = [];k = 1;answer = subreward(alpha,beta,p,gamma,n);armendfunction [answer] = subreward(alpha,beta,p,gamma,n)global arm kif n==1    [val,pull] = max([alpha/(alpha+beta),p]);    answer = val;    arm(k) = pull;    k=k+1;else    s = alpha/(alpha+beta) * (1 + gamma*subreward(alpha+1,beta,p,gamma,n-1));    f = beta/(alpha+beta) * (gamma*subreward(alpha,beta+1,p,gamma,n-1));    o = p*(1-gamma^n)/(1-gamma);    [val,pull] = max([s+f,o]);    answer = val;    arm(k) = pull;    k=k+1;endendHowever, this stores all 2^n - 1 values, rather than just n values.
Writing a program to find the optimal reward for a 2-armed Bernoulli bandit
dynamic programming