patchGNU Octave - Patches: patch #10016, [octave forge] (statistics)...

 
 

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

Submitter:  None
Submitted:  Sun 24 Jan 2021 10:10:24 AM UTC
   
 
Category:  Forge : new function Priority:  5 - Normal
Status:  Done Privacy:  Public
Assigned to:  philipnienhuis Originator Email:  -email is unavailable-
Open/Closed:  Closed
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 11 Feb 2021 09:10:35 PM UTC, comment #20: 

Thanks.
I added indents for the first two arguments, changed CRLF eOLs into just LF (also for patch #10019 BTW) and pushed the function + amended NEWS & Index.in under just your name here:
http://hg.code.sf.net/p/octave/statistics/rev/b641bfecec63

It'll be part of the next statistics package release just like slicesample.m). I have no idea hen that will happen, I'm not the maintainer, but given that quite a few new functions have been added and improvements have been made I have good hope that it won't be long.

Thank you very much for this contribution.

Closing report.

Philip Nienhuis <philipnienhuis>
Group Member
Wed 10 Feb 2021 10:56:51 PM UTC, comment #19: 

Several issues with the documentation. It contradicts itself, Inconsistent/undefined variable names and strings for the variable names were not provided. I have attached a new version.

(file #50852)

Anonymous
Wed 10 Feb 2021 08:49:20 AM UTC, comment #18: 

I got distracted with other occupations. In the mean time I fixed style issues (mostly spaces) and reorganized the Texinfo. All tests & demos pass and work, resp.
Could you please check if the attached is satisfactory for you?

(file #50848)

Philip Nienhuis <philipnienhuis>
Group Member
Mon 01 Feb 2021 04:30:54 AM UTC, comment #17: 

slicesample in see also should be acceptable now as I implemented it patch #10019

I don't think axis equal is better, the focus is on the sampling and the distribution, not the random function I chose to integrate. Axis equal makes the distribution very flat and does not make it clear the samples are mainly in the high probability area.

For the first question:

  • proprnd must always be given
  • One of pdf or logpdf must be given
  • One of proppdf or logproppdf must be given unless symmetric is true then neither need to be given
  • The remainder of the inputs are optional but will probably be given due to the way the sampler works.
Anonymous
Sun 31 Jan 2021 09:05:46 PM UTC, comment #16: 

Thanks.

Disclaimer:
I know next to nothing about mhsample, so all I can do is check the code for illogical, ambiguous and inefficient parts (found none), fix the texinfo header (usually the hardest part for new contributors; not that I'm much of a texinfo guru), and check style (found a few thingies). Plus check with on-line Matlab examples I can find but I suppose you've already done that (?)
Any functional bugs will have to be catched by users working with the function. (Perceive it this way: the Octave project doesn't really differ much there from many commercial SW vendors :-) )

While reviewing I got a few questions:

  1. Which property/value pairs are absolutely indispensable and which ones are optional?
  2. How about a statement 'axis equal' after the last mesh command in the first demo to make the plot look better?

As to seealso, I think we find it undesirable to refer to unimplemented functions there.

Philip Nienhuis <philipnienhuis>
Group Member
Sun 31 Jan 2021 02:18:55 AM UTC, comment #15: 

Replaced Heaviside with logical operators. I changed the output from no semicolon to fprintf on both demos. Plotted the function I'm integrating with the second demo.

The second see also, slicesample is not implemented. I think it should stay as a see also as that function does the same thing without needing to provide proppdf or proprnd. The paper Matlab references for slicesample has R code on the paper's author's website. https://www.cs.toronto.edu/~radford/homepage.html

The author suggests several different methods so getting the exact one Matlab uses may be difficult or is a method with the same inputs and outputs acceptable?

I'm the same person who submitted uniquetol. You can get my details from a patch submission in bug 59850.

(file #50811)

Anonymous
Sat 30 Jan 2021 08:04:41 PM UTC, comment #14: 

Nice work. Good adherence to GNU coding style, tests, demo.
Thanks very much!

I tried it out and found this:

Demo 1 gives a nice plot; nice!

-but:-

Demo 2 shows nothing ... no plot, no output. In addition it relies on the 'heaviside' function which is in the symbolic package. This implies adding a dependence on the symbolic package for the sake of just a demo. I don't think this is desirable.
Can you please come up with another 2nd demo? (or just drop it)

Once you have this sorted out I think I can push it.
Next then I'll probably push a few minor style & texinfo fixes - I'm sometimes a nitpicker, though surely not the most fanatical one out here :-)

Would you like to be credited for this contribution (well-deserved IMO)? In that case I can push it at least under your name, your email isn't required. We do need your name then.

(BTW your email can be handy for in case we get bug reports or other issues; but I don't know a secure way to keep it private.
It's easy to create an account here ("New user"), only a nickname and a valid email address is required. Visitors can't see it but unfortunately any registered Savannah user can.)

Philip Nienhuis <philipnienhuis>
Group Member
Sat 30 Jan 2021 02:46:35 PM UTC, comment #13: 

The complete mhsample function attached with documentation, tests, demos and following Octave style guide. I hope this is completely ready for inclusion into Octave.

(file #50809)

Anonymous
Fri 29 Jan 2021 11:58:45 AM UTC, comment #12: 

Looks like I've implemented mhsample completely with the second submission. I'll write the documentation and change the code style to line up with the style guide in the next few days. The code should not change from the second submission.

Anonymous
Fri 29 Jan 2021 11:16:05 AM UTC, comment #11: 

Uncommenting out L.4 works, script runs all the way to the end and I get a picture.
diary and pic attached.

I also tried with just L.6 commented out, but that gives:

:
[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 #50803,

Philip Nienhuis <philipnienhuis>
Group Member
Fri 29 Jan 2021 10:00:49 AM UTC, comment #10: 

Uncomment line 4. Comment out line 5. If that does not work do the same for line 6.

I have written mhsample for octave. I'm confident I have created the function correctly from the theory. As I do not have Matlab I can not be certain the format of the input and output are the same.

The layout of start and smpl are what I'm not sure about when the option nchain>1. What I need help with is what is the format of these variables.

Anonymous
Fri 29 Jan 2021 09:10:00 AM UTC, comment #9: 

You mean L.4, 5 or 6? which one do you want me to adapt?

Philip Nienhuis <philipnienhuis>
Group Member
Thu 28 Jan 2021 11:52:12 PM UTC, comment #8: 

Your test shows


start = rand (1, 2, nchain);


Is not a valid input. The Matlab documentation says "start is a row vector containing the start value of the Markov Chain, ..." which gives no indication of how different starting points should be input when nchain > 1 or if different starting points are possible in Matlab. I proposed two other potential inputs commented out with "Not sure how the starting point is input..." commented above them


start = rand (nchain, 2);

or

start = rand(1, 2);


Do either of the above work?

Once mhsample returns a result I need to compare the sizes of the output

size(smpl)
size(accept)


Again Matlab documentation is not clear about this for nchain>1. If I guessed correctly the code following mhsample should run fine and the values I calculate should line up with the answers, if not I expect an error or incorrect results.

The purpose of these tests is to work out what Matlab is doing that they have not documented.

Anonymous
Thu 28 Jan 2021 08:05:50 PM UTC, comment #7: 

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)

Philip Nienhuis <philipnienhuis>
Group Member
Tue 26 Jan 2021 11:57:51 AM UTC, comment #6: 

I decided to finish the m file implementation with the nchains propertie. I made guesses on how I think Matlab would take inputs and present the output. Unfortunately, Matlab's documentation does not give any examples using more than one chain. The original submission had no issues with all of Matlab's examples.

Metropolis-Hastings sampling is a method to generate random numbers from a probability distribution. It has the following advantages over other sampling methods:

  • The distribution does not need to be normalized.
  • It is particularly efficient in high dimensions where most of the pdf is near zero (no curse of dimensionality), selection rejection sampling will result in most values being rejected in this case

Its disadvantages are:

  • Points next to each other in a chain are correlated.
  • If the initial point is not in a location of high probability the initial samples are not useful.


I have included below an example of sampling from a distribution made from polynomials. I compare the sample mean, standard deviation along the axis direction (diagonal), and the covariance (off-diagonal) with what I calculated algebraically. This example demonstrated all of the inputs to the method and I have commented on what I am unsure about. In the file, there are also two demos. The first one samples a normal distribution and uses the sampled points to integrate the top half of a sphere. The second one samples a single variable normal distribution to find the normalization constant of a truncated normal distribution. I have tested all of Matlab's examples.


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);
close all;hold on;
[xx yy]=meshgrid(linspace(-1,1,25));
mesh(xx,yy,135/814*reshape(pdf([xx(:),yy(:)]),size(xx)),'facecolor','none');
plot(smpl(:,1,1),smpl(:,2,1),'.');
hold off;
mxy=mean(smpl,1);
sm=(smpl-mxy);
fprintf("Mean answer: [%f %f] \n",[-36/407 24/407]);
avg=mean(mxy,3)
fprintf("Diagonal of covariance answer: [%f %f] \n",[(1/3478629)*sqrt(1110349)*sqrt(3478629) (1/1159543)*sqrt(384653)*sqrt(1159543)]);
diagstd=mean(mean(sm.^2,1),3).^.5
fprintf("Off diagonal of covariance answer: %f \n",[-11346/165649]);
covar=mean(mean(prod(sm,2),1),3)
disp("Normally used to tune the method depends on proppdf and proprnd as well as pdf");
avgacc=mean(accept)%Normally used for tuning the method


If the implementation lines up with Matlab I will provide the help for the file and add the spaces to make it compliment with the Octave coding guidelines. As the c++ version did not significantly speed up the function (a few seconds faster for a 20-second runtime in pure Octave) I'm not sure if I should update the c++ version.

The new version is attached.

(file #50790)

Anonymous
Mon 25 Jan 2021 02:31:27 PM UTC, comment #5: 

Sure.
Please supply ready-baked instructions / code / data that I can feed 1:1 to Matlab. I'm unfamiliar with mhsample so can't figure that out myself.

Once mhsample works, please try the examples from the on-line Matlab docs. Doing that myself inspires most of my comments when reviewing external contributions :-)

Philip Nienhuis <philipnienhuis>
Group Member
Sun 24 Jan 2021 10:24:03 PM UTC, comment #4: 

It would be helpful if I could see the output with nchain set to 1, 2, 3 etc until a pattern could be observed. A small number of samples could be requested and the other parameters should no matter.

Anonymous
Sun 24 Jan 2021 10:00:07 PM UTC, comment #3: 

adapted title a bit

Philip Nienhuis <philipnienhuis>
Group Member
Sun 24 Jan 2021 09:58:37 PM UTC, comment #2: 

I have access to Matlab and Statistics toolbox, so if you have instructions for me I can try things out for you.

Philip Nienhuis <philipnienhuis>
Group Member
Sun 24 Jan 2021 11:28:51 AM UTC, comment #1: 

Thank you for the resubmission. 🙂

Kai Torben Ohlhus <siko1056>
Group Member
Sun 24 Jan 2021 10:10:24 AM UTC, original submission:  

In addition to what I said in bug #59924 I had not finished the implementation although all of the basic functionality was implemented.

I am/was hoping to get some information about the 'nchain' implementation, specifically, how is the output presented, is it an n-dimensional array or has Matlab interweaved the chains or something different. Also, does Matlab allow different starting positions for different chains, if so how is the input provided?

I'm also not sure what Matlab does for the default scaling of attempted step sizes when symmetric is true. I copied the example on the help page but I could imagine the step sizes being too big or small. This probably is not a big deal because proprnd provides this information if supplied.

Once I had finished the 'nchain' part I was then planning to write the documentation and change the code style to the Octave style guide (I find it harder to follow lines when there are many spaces so I like to add them after I have finished coding.)

I have been unable to post anything here from a little after the submission so I was unable to resubmit this earlier. I notice now it says a patch has been submitted. Have you finished implementing this function and written the documentation?

I have reattached the files

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #50852:  mhsample.m added by None (12KiB - text/plain)
file #50848:  mhsample.m added by philipnienhuis (12KiB - text/plain)
file #50811:  mhsample.m added by None (11KiB - text/plain)
file #50809:  mhsample.m added by None (10KiB - text/plain)
file #50803:  diary_L4.txt added by philipnienhuis (5KiB - text/plain)
file #50804:  diary_L4.png added by philipnienhuis (152KiB - image/png)
file #50790:  mhsample.m added by None (5KiB - text/plain)
file #50776:  mhsample.m added by None (4KiB - text/plain)
file #50777:  mhsample_ccc.cc added by None (2KiB - text/plain)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by philipnienhuis
  • -email is unavailable- added by siko1056 (Posted a comment)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only logged-in users can vote.

     

    Follow 16 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-02-11 philipnienhuis StatusIn Progress Done
        Open/ClosedOpen Closed
    2021-02-10 None Attached File- Added mhsample.m, #50852
    2021-02-10 philipnienhuis Attached File- Added mhsample.m, #50848
    2021-01-31 None Attached File- Added mhsample.m, #50811
    2021-01-30 philipnienhuis StatusNone In Progress
        Assigned toNone philipnienhuis
    2021-01-30 None Attached File- Added mhsample.m, #50809
    2021-01-29 philipnienhuis Attached File- Added diary_L4.txt, #50803
        Attached File- Added diary_L4.png, #50804
    2021-01-28 philipnienhuis Attached File- Added testcode_comment_6.m, #50801
    2021-01-26 None Attached File- Added mhsample.m, #50790
    2021-01-24 philipnienhuis Summary[octave forge] (statistics) mhsample implementation bug #59924 Resubmit as patch [octave forge] (statistics) mhsample implementation (over from bug #59924)
    2021-01-24 philipnienhuis Carbon-Copy- Added philipnienhuis
    2021-01-24 None Attached File- Added mhsample.m, #50776
        Attached File- Added mhsample_ccc.cc, #50777

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code