[Octave-patch-tracker] [patch #10016] [octave forge] (statistics) mhsample implementation (over from bug #59924)

[Philip Nienhuis]

Follow-up Comment #7, patch #10016 (project octave):

I adapted the code you entered in verbatim blocks a little, to make it
readable. You see, I'm no fan of Matlab code w/o any spaces and with a
compressed layout that makes it easier to follow the code. (Must be triggered
by former colleagues whose greater ambition it was to write complete regional
groundwater models in Matlab in just one line. Result: absolute illegible
gunk.)
Result is attached.

FYI Octave has a style guide, see here:
https://wiki.octave.org/Octave_style_guide

Anyway, running it in Matlab 2020b gives an error on L.20 (or L.19 in your
code):

>> diary on
>> nchain = 1000;                          % Number of chains
% Not sure how the starting point is input all three of the below work for my
implementation
% I allow different starting points, I hope Matlab does as well
%start=rand(nchain,2);
start = rand (1, 2, nchain);
>> %start=rand(1,2);%all chains start from same point
nsamples = 10000;                         % Length of chain
pdf = @(x) (-(x(:, 1) - 1 / 2) .^ 2 .* (x(:, 1) - 1) .* (x(:, 1) + 1) - ...
        x(:, 1) .* x(:, 2) - (x(:, 2) + 1/3) .^ 2 .* (x(:, 2) + 1) .* ...
        (x(:, 2) - 1) + 1) .* (x(:, 1) >= -1 & x(:, 1) <= 1 & ...
        x(:, 2) >= -1 & x(:, 2) <= 1);    % * 135 / 814; % The distribution we
wish to sample from, does not need to be normalized. Also assuming vectorized
function evaluation have different points down and different variables across
>> delta = .5;
proppdf = @(x, y) prod (unifpdf (y - x, -delta, delta), 2);  % pdf to choose
next point in chain
proprnd = @(x) x + rand (size (x)) * 2 * delta - delta;      %function to draw
random number from above pdf
>> sym = true;                               % if proppdf is a symmetric
distribution proppdf not used
K = 0;                                    % number of samples to discard at
the beginning
m = 2;                                    %disgard (m-1) every m samples
% not sure of the output sizes of smpl and accept. I have values of a chain
% dimension 1, variables dimension 2 and different chains dimension 3
[smpl, accept] = mhsample (start, nsamples, 'pdf', pdf, 'proppdf', proppdf,
...
                           'proprnd',proprnd, 'symmetric', sym, 'burnin', K,
...
                           'thin', m, 'nchain', nchain);
The logical indices in position 1 contain a true value outside of the array
bounds.

Error in mhsample (line 178)
    x0(acc,:) = y(acc,:); % preserves x's shape.
 
>> %% your original L.19 yields the same:
>>
[smpl,accept]=mhsample(start,nsamples,'pdf',pdf,'proppdf',proppdf,'proprnd',proprnd,'symmetric',sym,'burnin',K,'thin',m,'nchain',nchain);
The logical indices in position 1 contain a true value outside of the array
bounds.

Error in mhsample (line 178)
    x0(acc,:) = y(acc,:); % preserves x's shape.
 
>> 




(file #50801)
    _______________________________________________________

Additional Item Attachment:

File name: testcode_comment_6.m           Size:2 KB
    <https://file.savannah.gnu.org/file/testcode_comment_6.m?file_id=50801>



    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/patch/?10016>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/