How to create function handle using random output

Hi all,

I am trying to implement something called a "Tailored Randomized Block-MH" algorithm, which requires me to take a likelihood function, change a subset of parameters while fixing the others through the optimization step.

Suppose I have a function f([p,q,r,x,y,z]) and a randomization step tells me to hold [p,q,r] at [p_bar,q_bar,r_bar] fixed while changing [x,y,z]. I can type manually something like:

 maximand = @(x) f([p_bar,q_bar,r_bar,[x_1,x_2,x_3])

how would I do something like this using just the output from the output of, say, randperm(3)?

[Walter Roberson]

It is not clear to me what should be randomly permuted ?

[Jacob Thompson]

What is being randomly permuted is which variables are being changed and which are held fixed.

Suppose this function is f([p,q,r,x,y,z]), sometimes I have to modify p and z and hold q,r,x,y fixed, or possibly q,r,z get modified while p,x, and y are held fixed.

[Walter Roberson]

We need a bit more context. Is that random permutation being done on every call to the function? So for example,

fsolve(@(x) f(x(randperm(length(x))), x0)

which would randomize the order ever single time f was called?

Or are you wanting to do something like

function stuff
   determine randomly which variables to hold constant this time
   fun = @(x) f(mix of held and non-held parameters using x)
   fsolve(fun, x0)
end

where within any one call to the function such a fsolve() or fmincon() the called function acts the same way, but you are in some sort of situation where that optimizer is being called repeatedly ?

[Jacob Thompson]

function stuff
determine randomly which variables to hold constant this time
fun = @(x) f(mix of held and non-held parameters using x)
fsolve(fun, x0)
end

This seems much closer to what I'm looking for, yes. And yes for the algorithm I am trying to implement, an optimizer is called repeatedly, and each time it is called I want to vary which parameters are fixed at the previous optima and which are allowed to vary.

[Walter Roberson]

Which specific optimizer are you using? Some of them make it easier than others.

[Jacob Thompson]

I usually use fminunc but anything capable of simulated annealing will work

Answer 1

[Jan]

Maybe:

 maximand = @(x) f([p_bar, q_bar, r_bar, x(randperm(3, 3))])

[Walter Roberson]

That potentially returns one of six different values, depending on which permutation occurs. It might possibly make sense to me to deliberately try all permutations of x, six calls to f, but I am having difficulty thinking of a context in which randomly selecting would be useful.

[Jacob Thompson]

That's not quite what I'm looking for. I don't want the order of variables to be randomized once we've already decided which parameters are fixed and which are allowed to vary to maximize the maximand, but rather I want to randomize the choice of which parameters to fix and which to vary in order to maximize the maximand.

Answer 2

[Walter Roberson]

For fmincon and simulannealbnd the easiest way to handle this is to use the same function call in each case, but set the ub and lb the same for the entries that are to be fixed for this run.

For example,

rvars = sort(randperm(numel(UB),NumFixedVars));

fixedvals = LB(rvars) + rand(1,NumFixedVars) .* (UB(rvars) - LB(rvars));

LB(rvars) = fixedvals; UB(rvars) = fixedvals;
x0(rvars) = fixedvals;
[bestx, fvals, exitflag] = fmincon(fun, x0, A, b, Aeq, beq, LB, UB, nonlcon, opts);

This code picks random variable indices. Then it picks random values between the lower and upper bound for those variables, and sets the lower and upper bound to be the same for those variables, and proceeds to run the optimization.

You might, of course, have had completely different logic in mind as to how to choose the fixed values for the variables, so alter the fixedvals assignment as appropriate for your situation.

Answer 3

[Jacob Thompson]

I'm quite certain that I figured this out, thanks for the help guys.

function liki = fmin(x,current,block)
% y = zeros(1,16);
noblock = setdiff(1:16,block);
y = zeros(size(current));
y(noblock) = current(noblock);
elements = length(block);
for i = 1:elements
    y(block(i)) = x(i);
end
liki = -dsgeliki(y) - prior(y);
end

[Walter Roberson]

Vectorized:

 y(block) = x(1:length(block)) ;