{"id":2351,"date":"2010-03-01T07:03:41","date_gmt":"2010-03-01T06:03:41","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=2351"},"modified":"2010-09-23T12:16:12","modified_gmt":"2010-09-23T11:16:12","slug":"an-alternative-to-binopdf-using-the-nag-toolbox-for-matlab","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=2351","title":{"rendered":"An alternative to binopdf using the NAG toolbox for MATLAB"},"content":{"rendered":"<h3><strong><strong>The Problem<br \/>\n<\/strong><\/strong><\/h3>\n<p>A MATLAB user came to me with a complaint recently.\u00a0 He had a piece of code that made use of the <a href=\"http:\/\/www.mathworks.co.uk\/products\/statistics\/\">MATLAB Statistics toolbox<\/a> but couldn&#8217;t get reliable access to\u00a0 a stats toolbox license from my employer&#8217;s server.\u00a0 You see, although we have hundreds of licenses for MATLAB itself, we only have a few dozen licenses for the statistics toolbox.\u00a0 This has always served us well in the past but use of this particular toolbox is on the rise and so we sometimes run out which means that users are more likely to get the following error message these days<\/p>\n<pre>??? License checkout failed.\r\nLicense Manager Error -4\r\nMaximum number of users for Statistics_Toolbox reached.\r\nTry again later.\r\nTo see a list of current users use the lmstat utility or contact your License Administrator.<\/pre>\n<p>The user had some options if he wanted to run his code via our shared network licenses for MATLAB (rather than rewrite in a free language or buy his own, personal license for MATLAB):<\/p>\n<ul>\n<li>Wait until a license became free.<\/li>\n<li>Wait until we found funds for more licenses.<\/li>\n<li>Buy an additional license token for himself and add it to the network (235 pounds +VAT currently)<\/li>\n<li>Allow me to gave a look at his code and come up with another way<\/li>\n<\/ul>\n<p>He went for the final option.\u00a0 So, I sat with a colleague of mine and we looked through the code.\u00a0 It turned out that there was only one line in the entire program that used the Statistics toolbox and all that line did was call binopdf.\u00a0 That one function call almost cost him 235 quid!<\/p>\n<p>Now, <a href=\"http:\/\/www.mathworks.com\/access\/helpdesk\/help\/toolbox\/stats\/binopdf.html\">binopdf()<\/a> simply evaluates the <a href=\"http:\/\/en.wikipedia.org\/wiki\/Binomial_distribution\">binomial probability density function<\/a> and the definition doesn&#8217;t look too difficult so you may wonder why I didn&#8217;t just cook up my own version of binopdf for him?\u00a0 Quite simply, I don&#8217;t trust myself to do as good a job as The Mathworks.\u00a0 There&#8217;s bound to be gotchas and I am bound to not know them at first.\u00a0 What would you rather use, binpodf.m from Mathworks or binopdf.m from &#8216;Dodgy Sysadmin&#8217;s Functions Inc.&#8217;?\u00a0 Your research depends upon this function remember&#8230;<\/p>\n<h3><strong><strong>NAG to the rescue!<br \/>\n<\/strong><\/strong><\/h3>\n<p>I wanted to use a version of this function that was written by someone I trust so my mind immediately turned to NAG (<a href=\"http:\/\/www.nag.co.uk\/\">Numerical Algorithms Group<\/a>) and their <a href=\"http:\/\/www.nag.co.uk\/numeric\/MB\/start.asp\">MATLAB toolbox<\/a> which my employer has a full site license for.\u00a0 The NAG Toolbox for MATLAB has a function called <strong>g01bj<\/strong> which does what we need and a whole lot more.\u00a0 The basic syntax is<\/p>\n<pre>[plek, pgtk, peqk, ifail] = g01bj(n, p, x)<\/pre>\n<p><strong> <\/strong><\/p>\n<p>From the NAG documentation: Let X denote a random variable having a binomial distribution with parameters n and p.  <strong>g01bj <\/strong>calculates the following probabilities<\/p>\n<ul>\n<li>plek = Prob{X&lt;=x}<\/li>\n<li>pgtk = Prob{X &gt; x}<\/li>\n<li>peqk = Prob{X = x}<\/li>\n<\/ul>\n<p>MATLAB&#8217;s binopodf only calculates Prob(X=x) so the practical upshot is that for the following MATLAB code<\/p>\n<pre>n=200;\r\nx=0;\r\np=0.02;\r\ny=binopdf(x,n,p)<\/pre>\n<p>the corresponding NAG Toolbox code (on a 32bit machine) is<\/p>\n<pre>n=int32(200);\r\nx=int32(0);\r\np=0.02;\r\n[plek, pgtk, y, ifail] = g01bj(n, p, x);<\/pre>\n<p>both of these code snippets gives<\/p>\n<pre>y =\r\n    0.0176<\/pre>\n<p>Further tests show that NAG&#8217;s and MATLAB&#8217;s results agree either exactly or to within what is essentially numerical noise for a range of values.  The only thing that remained was to make NAG&#8217;s function look more &#8216;normal&#8217; for the user in question which meant creating a file called nag_binopdf.m which contained the following<\/p>\n<pre>function pbin = nag_binopdf(x,n,p)\r\n   x=int32(x);\r\n   n=int32(n);\r\n   [~,~,pbin,~] = g01bj(n,p,x);\r\nend<\/pre>\n<p>Now the user had a function called nag_binopdf that behaved just like the Statistics toolbox&#8217;s binopdf, <strong>for his particular application<\/strong>, argument order and everything.<\/p>\n<h3><strong>It&#8217;s not all plain sailing though!<\/strong><\/h3>\n<p>As good as the NAG toolbox is, it&#8217;s not perfect.\u00a0 Note that I placed particular emphasis on the fact that we had come up with a suitable, drop-in replacement for binopdf for this <strong>user&#8217;s particular application<\/strong>.\u00a0 The salient points of his application are:<\/p>\n<ul>\n<li><strong>His input argument, x, is just a single value &#8211; not a vector of values.<\/strong> The NAG library wasn&#8217;t written with MATLAB in mind, it&#8217;s written in Fortran, and so many of its functions are not vectorised (some of them are though).\u00a0 In Fortran this is no problem at all but in MATLAB it&#8217;s a pain.\u00a0 I could, of course, write a nag_binopodf that <strong>looks<\/strong> like it&#8217;s vectorised by putting g01bj in a loop hidden away from the user but I&#8217;d soon get found out because the performance would suck.\u00a0 To better support MATLAB users, NAG needs to write properly vectorised versions of its functions.<\/li>\n<li><strong>He only wanted to run his code on 32bit machines.<\/strong> If he had wanted to run it on a 64bit machine then he would have to change all of those <strong>int32() <\/strong>calls to <strong>int64()<\/strong>.\u00a0 I agree that this is no big deal but it&#8217;s something that he wouldn&#8217;t have to worry about if he had coughed up the cash for a Mathworks statistics license.<strong> <\/strong><\/li>\n<\/ul>\n<h3>Final thoughts<\/h3>\n<p>I use the NAG toolbox for MATLAB to do this sort of thing all the time and have saved researchers at my university several thousand pounds over the last year or so through unneeded Mathworks toolbox licenses.\u00a0 As well as statistics, the NAG toolbox is also good for optimisation (local and global), curve fitting, splines, wavelets, finance,partial differential equations and probably more.\u00a0 Not everyone likes this approach and it is not always suitable but when it works, it works well.\u00a0 Sometimes, you even get a significant speed boost into the bargain (check out my <a href=\"..\/?p=1552\">interp1 article<\/a> for example).<\/p>\n<p>Comments are, as always, welcomed.<\/p>\n<p><strong>Blog posts similar to this one<\/strong><\/p>\n<ul>\n<li><a href=\"https:\/\/www.walkingrandomly.com\/?p=1488\">A faster version of MATLAB&#8217;s fsolve using the NAG Toolbox for MATLAB<\/a><\/li>\n<li><a href=\"https:\/\/www.walkingrandomly.com\/?p=2782\">An alternative to the ranksum function using the NAG Toolbox for MATLAB<\/a><\/li>\n<li><a href=\"https:\/\/www.walkingrandomly.com\/?p=2678\">A free version of the pdist command for MATLAB<\/a><\/li>\n<\/ul>\n","protected":false},"excerpt":{"rendered":"<p>The Problem A MATLAB user came to me with a complaint recently.\u00a0 He had a piece of code that made use of the MATLAB Statistics toolbox but couldn&#8217;t get reliable access to\u00a0 a stats toolbox license from my employer&#8217;s server.\u00a0 You see, although we have hundreds of licenses for MATLAB itself, we only have a [&hellip;]<\/p>\n","protected":false},"author":3,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"footnotes":"","jetpack_publicize_message":"","jetpack_publicize_feature_enabled":true,"jetpack_social_post_already_shared":false,"jetpack_social_options":{"image_generator_settings":{"template":"highway","default_image_id":0,"font":"","enabled":false},"version":2}},"categories":[4,11,28,7],"tags":[],"class_list":["post-2351","post","type-post","status-publish","format-standard","hentry","category-math-software","category-matlab","category-nag-library","category-programming"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p3swhs-BV","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/2351","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=2351"}],"version-history":[{"count":16,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/2351\/revisions"}],"predecessor-version":[{"id":2369,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/2351\/revisions\/2369"}],"wp:attachment":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=2351"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=2351"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=2351"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}