turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Find a Community

- Home
- /
- Analytics
- /
- Stat Procs
- /
- Matrix functions in Proc MCMC

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

07-05-2015 04:44 PM

Hello everyone,

I have noticed that PROC MCMC is limited to apply functions on matrices. We can add, subtract, multiply matrices, or take determinant and interse of them. But if I want exponential or log of each entry of a matrix, I need to program it. For example:

array MAT[1000,1000] ;

array cov[1000,1000] ;

do i = 1 to 1000;

do j = i+1 to 1000;

cov[i, j] = exp( MAT[i, j]);

cov[j, i] = cov[i, j];

end;

cov[i,i] = exp(MAT[i, i]);

end;

This slows down the program. Don't you know any simpler or more efficient way of applying functions on matrices? How can I take the exponential of a matrix without a loop?

Accepted Solutions

Solution

07-06-2015
05:38 PM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Mehdi_R

07-06-2015 05:38 PM

You can make use of the matrixfunctions in fcmp. But unfortunately they can not be called immediately from a datastep or procedure. Therefore, what you can do is to make a matrix-function which call the fcmp's expmatrix-function. Your function can then be called from a datastep or from a procedure.

proc fcmp outlib = work.func.matrix;

subroutine matrixexp(m[*,*],t,y[*,*]);

outargs y;

call expmatrix(m, t, y);

endsub;

run;

option cmplib=(work.func);

data _NULL_;

array a{2,2} _temporary_;

array y{2,2} _temporary_;

call matrixexp(a,4,y);

run;

I haven't tested if the functions declared in this way can be called from PROC MCMP, but they work within some other procedures (nlmixed for instance).

I made an idea on this site about making the matrix functions directly available without first declare them with PROC FCMP. To my surprise it didnt get so positive votes :smileyshocked: https://communities.sas.com/ideas/1570

All Replies

Solution

07-06-2015
05:38 PM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Mehdi_R

07-06-2015 05:38 PM

You can make use of the matrixfunctions in fcmp. But unfortunately they can not be called immediately from a datastep or procedure. Therefore, what you can do is to make a matrix-function which call the fcmp's expmatrix-function. Your function can then be called from a datastep or from a procedure.

proc fcmp outlib = work.func.matrix;

subroutine matrixexp(m[*,*],t,y[*,*]);

outargs y;

call expmatrix(m, t, y);

endsub;

run;

option cmplib=(work.func);

data _NULL_;

array a{2,2} _temporary_;

array y{2,2} _temporary_;

call matrixexp(a,4,y);

run;

I haven't tested if the functions declared in this way can be called from PROC MCMP, but they work within some other procedures (nlmixed for instance).

I made an idea on this site about making the matrix functions directly available without first declare them with PROC FCMP. To my surprise it didnt get so positive votes :smileyshocked: https://communities.sas.com/ideas/1570

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to JacobSimonsen

07-07-2015 02:57 AM

Thanks Jacob, but I think still there is problem. Please take a look at this page:

Base SAS(R) 9.2 Procedures Guide

Good news is that it works in PROC MCMC. But bad news is that the page above says that CALL EXPMATRIX does not exponentiate each element of the matrix! While I need to exponentiate each element of a matrix in PROC MCMC. Do you have another suggestion or comment?

Thanks

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Mehdi_R

07-07-2015 02:15 AM

It was awesome! Thanks.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Mehdi_R

07-07-2015 03:19 AM

Sorry, I did not read your question correct. If you want to exponentiate elementwise, then change the fcmp function to

proc fcmp outlib = work.func.matrix;

subroutine matrixexp(m[*,*],y[*,*]);

outargs y;

do i=1 to dim(m,1);

do j=1 to dim(m,2);

y[i,j]=exp(m[i,j]);

end;

end;

endsub;

run;

Btw, you can not avoid the loop. Therefore hiding the matrix exponential within a function will not speed up the program.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to JacobSimonsen

07-07-2015 04:12 AM

Thanks Jacob for your quick reply.

I thought that there might be a more efficient way to calculate exponential of a matrix (matrix of exponential of entries) in PROC MCMC. Because MCMC involves heavy calculations and I wanted to avoid such loop... It seems there is no built-in function in SAS which calculates exponential of a matrix...

Do you think writing loops in proc mcmc is faster or writing loops in proc fcmp and then call the subroutine to proc mcmc?

Thanks again Jacob

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Mehdi_R

07-07-2015 04:37 AM

No, I think it will be almost same speed. Also, even if there was a built in function to make elementwise exponential, then there still is the loop;

Though, it could be made such that the elements are exponentiated in parallel, but don't ask me how that should be programmed

If you want to gain a bit, you can add m to outargs. I think some copying of values then can be avoided. It will not change much.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to JacobSimonsen

07-07-2015 04:47 AM

Thank you Jacob. It was really helpful