In a recent article, Citing software in research papers, I discussed how to cite various software packages. One of the commentators suggested that I should contact The Mathworks if I wanted to know how to cite MATLAB. I did this and asked for permission to blog the result. This is what they suggested:
To cite MATLAB (and in this case a toolbox) you can use this:
MATLAB and Statistics Toolbox Release 2012b, The MathWorks, Inc., Natick, Massachusetts, United States.
This is the format prescribed by both the Chicago Manual of Style and the Microsoft Manual of Style. If you use a different style guide, you may apply a different format, but should observe the capitalization shown above, and include the appropriate release number. If the MathWorks release number (in the format YYYYa or YYYYb) is not readily available, you can use the point release numbers for the software, as in:
MATLAB 8.0 and Statistics Toolbox 8.1, The MathWorks, Inc., Natick, Massachusetts, United States.
Thanks to The Mathworks for allowing me to reproduce this communication here at WalkingRandomly.
xkcd is a popular webcomic that sometimes includes hand drawn graphs in a distinctive style. Here’s a typical example
In a recent Mathematica StackExchange question, someone asked how such graphs could be automatically produced in Mathematica and code was quickly whipped up by the community. Since then, various individuals and communities have developed code to do the same thing in a range of languages. Here’s the list of examples I’ve found so far
- xkcd style graphs in Mathematica. There is also a Wolfram blog post on this subject.
- xkcd style graphs in R. A follow up blog post at vistat.
- xkcd style graphs in LaTeX
- xkcd style graphs in Python using matplotlib
- xkcd style graphs in MATLAB. There is now some code on the File Exchange that does this with your own plots.
- xkcd style graphs in Euler
Any I’ve missed?
Back when Mathematica 8 was released I tried to work out how many MATLAB toolboxes you’d need to buy to have the same functionality and came up with 9 toolboxes. Readers of WalkingRandomly suggested several more in the comments. Now that Mathematica 9 has been released, I thought I’d work through the exercise again.
So I think that Mathematica 9 contains at least some of the functionality of the following 18 MATLAB toolboxes. Click on the relevant toolbox for more information or an example.
- Control systems toolbox
- Curve fitting toolbox
- Datafeed toolbox
- Database toolbox
- DSP System toolbox
- Econometrics toolbox
- Finance toolbox
- Financial Instruments toolbox
- Global Optimization Toolbox
- Image Acquisition Toolbox
- Image Processing toolbox
- MATLAB Compiler
- Optimization toolbox
- Parallel Computing toolbox
- Signal Processing toolbox
- Statistics toolbox
- Symbolic toolbox
- Wavelet toolbox
I use both Mathematica and MATLAB extensively and sincerely wish that MATLAB had this level of integration. Does anyone have evidence of any I might have missed (or shouldn’t have included)?
MATLAB Mobile has been around for Apple devices for a while now but Android users have had to make do with third party alternatives such as MATLAB Commander and MLConnect. All that has now changed with the release of MATLAB Mobile for Android.
MATLAB Mobile is NOT MATLAB running on your phone
While MATLAB Mobile is a very nice and interesting product, there is one thing you should get clear in your mind– this is not a full version of MATLAB on your phone or Tablet. MATLAB Mobile is essentially a thin client that connects to an instance of MATLAB running on your desktop or The Mathworks Cloud. In other words, it doesn’t work at all if you don’t have a network connection or a licensed copy of MATLAB.
What if you do want to run MATLAB code directly on your phone?
While it is unlikely that we’ll see a full version of MATLAB compiled for Android devices any time soon, Android toting MATLABers have a couple of other options available to them in addition to MATLAB Mobile.
- Octave for Android Octave is a free, open source alternative to MATLAB that can run many .m file scripts and functions. Corbin Champion has ported it to Android and although it is still a work in progress, it works very well.
- Mathmatiz – Small and light, this free app understands a subset of the MATLAB language and can do basic plotting.
- Addi – Much smaller and less capable than Octave for Android, this is Corbin Champion’s earlier attempt at bringing a free MATLAB clone to Android. It is based on the Java library, JMathLib.
Simulink from The Mathworks is widely used in various disciplines. I was recently asked to come up with a list of alternative products, both free and commercial.
Here are some alternatives that I know of:
- MapleSim – A commercial Simuink replacement from the makers of the computer algebra system, Maple
- OpenModelica -An open-source Modelica-based modeling and simulation environment intended for industrial and academic usage
- Wolfram SystemModeler – Very new commercial product from the makers of Mathematica. Click here for Wolfram’s take on why their product is the best.
- xcos – This free Simulink alternative comes with Scilab.
I plan to keep this list updated and, eventually, include more details. Comments, suggestions and links to comparison articles are very welcome. If you have taught a course using one of these alternatives and have experiences to share, please let me know. Similarly for anyone who was switched (or attempted to switch) their research from Simulink. Either comment to this post or contact me directly.
I’ve nothing against Simulink but would like to get a handle on what else is out there.
I felt like playing with Julia and MATLAB this Sunday morning. I found some code that prices European Options in MATLAB using Monte Carlo simulations over at computeraidedfinance.com and thought that I’d port this over to Julia. Here’s the original MATLAB code
function V = bench_CPU_European(numPaths) %Simple European steps = 250; r = (0.05); sigma = (0.4); T = (1); dt = T/(steps); K = (100); S = 100 * ones(numPaths,1); for i=1:steps rnd = randn(numPaths,1); S = S .* exp((r-0.5*sigma.^2)*dt + sigma*sqrt(dt)*rnd); end V = mean( exp(-r*T)*max(K-S,0) )
I ran this a couple of times to see what results I should be getting and how long it would take for 1 million paths:
tic;bench_CPU_European(1000000);toc V = 13.1596 Elapsed time is 6.035635 seconds. >> tic;bench_CPU_European(1000000);toc V = 13.1258 Elapsed time is 5.924104 seconds. >> tic;bench_CPU_European(1000000);toc V = 13.1479 Elapsed time is 5.936475 seconds.
The result varies because this is a stochastic process but we can see that it should be around 13.1 or so and takes around 6 seconds on my laptop. Since it’s Sunday morning, I am feeling lazy and have no intention of considering if this code is optimal or not right now. I’m just going to copy and paste it into a julia file and hack at the syntax until it becomes valid Julia code. The following seems to work
function bench_CPU_European(numPaths) steps = 250 r = 0.05 sigma = .4; T = 1; dt = T/(steps) K = 100; S = 100 * ones(numPaths,1); for i=1:steps rnd = randn(numPaths,1) S = S .* exp((r-0.5*sigma.^2)*dt + sigma*sqrt(dt)*rnd) end V = mean( exp(-r*T)*max(K-S,0) ) end
I ran this on Julia and got the following
julia> tic();bench_CPU_European(1000000);toc() elapsed time: 36.259000062942505 seconds 36.259000062942505 julia> bench_CPU_European(1000000) 13.114855104505445
The Julia code appears to be valid, it gives the correct result of 13.1 ish but at 36.25 seconds is around 6 times slower than the MATLAB version. The dog needs walking so I’m going to think about this another time but comments are welcome.
Update (9pm 7th October 2012): I’ve just tried this Julia code on the Linux partition of the same laptop and 1 million paths took 14 seconds or so:
tic();bench_CPU_European(1000000);toc() elapsed time: 14.146281957626343 seconds
I built this version of Julia from source and so it’s at the current bleeding edge (version 0.0.0+98589672.r65a1 Commit 65a1f3dedc (2012-10-07 06:40:18). The code is still slower than the MATLAB version but better than the older Windows build
Update: 13th October 2012
Over on the Julia mailing list, someone posted a faster version of this simulation in Julia
function bench_eu(numPaths) steps = 250 r = 0.05 sigma = .4; T = 1; dt = T/(steps) K = 100; S = 100 * ones(numPaths,1); t1 = (r-0.5*sigma.^2)*dt t2 = sigma*sqrt(dt) for i=1:steps for j=1:numPaths S[j] .*= exp(t1 + t2*randn()) end end V = mean( exp(-r*T)*max(K-S,0) ) end
On the Linux partition of my test machine, this got through 1000000 paths in 8.53 seconds, very close to the MATLAB speed:
julia> tic();bench_eu(1000000);toc() elapsed time: 8.534484148025513 seconds
It seems that, when using Julia, one needs to unlearn everything you’ve ever learned about vectorisation in MATLAB.
Update: 28th October 2012
Members of the Julia team have been improving the performance of the randn() function used in the above code (see here and here for details). Using the de-vectorised code above, execution time for 1 million paths in Julia is now down to 7.2 seconds on my machine on Linux. Still slower than the MATLAB 2012a implementation but it’s getting there. This was using Julia version 0.0.0+100403134.r0999 Commit 099936aec6 (2012-10-28 05:24:40)
tic();bench_eu(1000000);toc() elapsed time: 7.223690032958984 seconds 7.223690032958984
- Laptop model: Dell XPS L702X
- CPU:Intel Core i7-2630QM @2Ghz software overclockable to 2.9Ghz. 4 physical cores but total 8 virtual cores due to Hyperthreading.
- RAM: 8 Gb
- OS: Windows 7 Home Premium 64 bit and Ubuntu 12.04
- MATLAB: 2012a
- Julia: Original windows version was Version 0.0.0+94063912.r17f5, Commit 17f50ea4e0 (2012-08-15 22:30:58). Several versions used on Linux since, see text for details.
Pop quiz: What does the following line of MATLAB code do?
If you said ‘It changes the seed of the random number generator to 10′ you get half a point.
‘Only half a point!?’ I hear you say accusingly ‘but it says so in my book [for example, 1-3], why not a full point?’
You only get a full point if you’d said something like ‘It changes the seed of the random number generator to 10 and it also changes the random number generator from the high quality, default Mersenne Twister generator to a lower quality legacy random number generator.‘
OK, how about this one?
This behaves in a very similar manner– it changes both the seed and the type of the underlying generator. However, the random number generator it switches to this time is an even older one that was introduced as far back as MATLAB version 4. It is not very good at all by modern standards!
A closer look
Open up a fresh copy of a recent version of MATLAB and ask it about the random number generator it’s using
>> RandStream.getGlobalStream ans = mt19937ar random stream (current global stream) Seed: 0 NormalTransform: Ziggurat
mt1993ar refers to a particular variant of the Mersenne Twister algorithm– an industry strength random number generator that’s used in many software packages and simulations. It’s been the default generator in MATLAB since 2007a. Change the seed using the modern (since 2011a), recommended syntax and ask again:
>> rng(10) >> RandStream.getGlobalStream ans = mt19937ar random stream (current global stream) Seed: 10 NormalTransform: Ziggurat
This is behaving exactly as you’d expect, you ask it to change the seed and it changes the seed…nothing more, nothing less. Now, let’s use the older syntax
>> rand('state',10) >> RandStream.getGlobalStream ans = legacy random stream (current global stream) RAND algorithm: V5 (Subtract-with-Borrow), RANDN algorithm: V5 (Ziggurat)
The random number generator has completely changed! We are no longer using the Mersenne Twister algorithm, we are now using a ‘subtract with borrow’ [see reference 4 for implementation details] generator which has been shown to have several undesirable issues [5-7].
Let’s do it again but this time using the even older ‘seed’ version:
>> rand('seed',10) >> RandStream.getGlobalStream ans = legacy random stream (current global stream) RAND algorithm: V4 (Congruential), RANDN algorithm: V5 (Ziggurat)
Now, this random number generator is ancient by computing standards. It also has a relatively tiny period of only 2 billion or so. For details see 
Why this matters
Now, all of this is well documented so you may wonder why I am making such a big issue out of it. Here are my reasons
- I often get sent MATLAB code for the purposes of code-review and optimisation. I see the old seeding syntax a LOT and the program’s authors are often blissfully unaware of the consequnces.
- The old syntax looks like all it should do is change the seed. It doesn’t! Before 2007a, however, it did!
- The old syntax is written in dozens of books because it was once the default, correct syntax to use.
- Many users don’t read the relevent section of the MATLAB documentation because they have no idea that there is a potential issue. They read a book or tutorial..it says to use rand(‘state’,10) so they do.
- MATLAB doesn’t use the old generators by default any more because they are not very good [4-7]!
- Using these old generators may adversely affect the quality of your simulation.
The bottom line
Don’t do either of these to change the seed of the default generator to 10:
Do this instead:
Only if you completely understand and accept the consequences of the older syntax should you use it.
1. ‘MATLAB – A practical introduction to programming and problem solving’, 2009,Stormy Attaway
2. MATLAB Guide (Second Edition), 2005, Desmond Higham and Nicholas Higham
3. Essential MATLAB for Engineers and Scientists (Fourth Edition), 2009, Hahn and Valentine
4. Numerical Computing with MATLAB, 2004, Cleve Moler (available online)
5. Why does the random number generator in MATLAB fail a particular test of randomness? The Mathworks, retreived 26th September 2012
6. A strong nonrandom pattern in Matlab default random number generator, 2006, Petr Savicky, retreived 26th September 2012
7. Learning Random Numbers: A Matlab Anomaly, 2008, Petr Savicky and Marko Robnik-Šikonja, Applied Artificial Intelligence, Vol22 issue 3, pp 254-265
Other posts on random numbers in MATLAB
>> A=[2 1 3;-1 0 7; 0 -1 -1]; >> [Q R]=qr(A) Q = -0.8944 -0.1826 0.4082 0.4472 -0.3651 0.8165 0 0.9129 0.4082 R = -2.2361 -0.8944 0.4472 0 -1.0954 -4.0166 0 0 6.5320
Variable precision matrices do not (I’m using MATLAB 2012a and the symbolic toolbox here).
>> a=vpa([2 1 3;-1 0 7; 0 -1 -1]); >> [Q R]=qr(a) Undefined function 'qr' for input arguments of type 'sym'.
It turns out that MATLAB and the symbolic toolbox CAN do variable precision QR factorisation….it’s just hidden a bit. The following very simple function, vpa_qr.m, shows how to get at it
function [ q,r ] = vpa_qr( x ) result = feval(symengine,'linalg::factorQR',x); q=result(1); r=result(2); end
Let’s see how that does
>> a=vpa([2 1 3;-1 0 7; 0 -1 -1]); >> [Q R]=vpa_qr(a);
I’ve suppressed the output because it’s so large but it has definitely worked. Let’s take a look at the first element of Q for example
>> Q(1) ans = 0.89442719099991587856366946749251
Which is correct to the default number of variable precision digits, 32. Of course we could change this to anything we like using the digits function.
- Download vpa_qr.m (needs symbolic toolbox)
The Mathworks sell dozens of toolboxes for MATLAB but they are not the only ones doing it. Many 3rd party vendors also sell MATLAB toolboxes and many of them are of extremely high quality. Here are some of my favourites
- The NAG Toolbox for MATLAB – Over 1500 high quality numerical routines from the Numerical Algorithms Group. Contains local and global optimisation, statistics, finance, wavelets, curve fitting, special functions and much much more. Superb!
- AccelerEyes Jacket – Very fast and comprehensive GPU computing toolbox for MATLAB. I’ve used it a lot.
- Multi-precision toolbox for MATLAB – A colleague at Manchester recently told me about this one as his group uses it a lot in their research. I’ve yet to play with it but it looks great.
Which commercial, 3rd party toolboxes do you use/rate and why?
If you like this list, you may also like my list of high quality, free MATLAB toolboxes
A MATLAB user at Manchester University contacted me recently asking about Black-Scholes option pricing. The MATLAB Financial Toolbox has a range of functions that can calculate Black-Scholes put and call option prices along with several of the sensitivities (or ‘greeks‘) such as blsprice, blsdelta and so on.
The user’s problem is that we don’t have any site-wide licenses for the Financial Toolbox. We do, however, have a full site license for the NAG Toolbox for MATLAB which has a nice set of option pricing routines. Even though they calculate the same things, NAG Toolbox option pricing functions look very different to the Financial Toolbox ones and so I felt that a Rosetta Stone type article might be useful.
For Black-Scholes option pricing, there are three main differences between the two systems:
- The Financial Toolbox has separate functions for calculating the option price and each greek (e.g. blsprice, blsgamma, blsdelta etc) whereas NAG calculates the price and all greeks simultaneously with a single function call.
- Where appropriate, The MATLAB functions calculate Put and Call values with one function call whereas with NAG you need to explicitly specify Call or Put.
- NAG calculates more greeks than MATLAB.
The following code example pretty much says it all. Any variable calculated with the NAG Toolbox is prefixed NAG_ whereas anything calculated with the financial toolbox is prefixed MW_. When I developed this, I was using MATLAB 2012a with NAG Toolbox Mark 22.
%Input parameters for both NAG and MATLAB. Price=50; Strike=40; Rate=0.1; Time=0.25; Volatility=0.3; Yield=0; %calculate all greeks for a put using NAG [NAG_Put, NAG_PutDelta, NAG_Gamma, NAG_Vega, NAG_PutTheta, NAG_PutRho, NAG_PutCrho, NAG_PutVanna,... NAG_PutCharm, NAG_PutSpeed, NAG_PutColour, NAG_PutZomma,NAG_PutVomma, ifail] =... s30ab('p', Strike, Price, Time, Volatility, Rate, Yield); %calculate all greeks for a Call using NAG [NAG_Call, NAG_CallDelta, NAG_Gamma, NAG_Vega, NAG_CallTheta, NAG_CallRho, NAG_CallCrho, NAG_CallVanna,... NAG_CallCharm, NAG_CallSpeed, NAG_CallColour,NAG_CallZomma, NAG_CallVomma, ifail] = ... s30ab('c', Strike, Price, Time, Volatility, Rate, Yield); %Calculate the Elasticity (Lambda) NAG_CallLambda = Price/NAG_Call*NAG_CallDelta; NAG_PutLambda = Price/NAG_Put*NAG_PutDelta; %Calculate the same set of prices and greeks using the MATLAB Finance Toolbox [MW_Call, MW_Put] = blsprice(Price, Strike, Rate, Time, Volatility, Yield); [MW_CallDelta, MW_PutDelta] = blsdelta(Price, Strike, Rate, Time, Volatility, Yield); MW_Gamma = blsgamma(Price, Strike, Rate, Time, Volatility, Yield); MW_Vega = blsvega(Price, Strike, Rate, Time, Volatility, Yield); [MW_CallTheta, MW_PutTheta] = blstheta(Price, Strike, Rate, Time,Volatility, Yield); [MW_CallRho, MW_PutRho]= blsrho(Price, Strike, Rate, Time, Volatility,Yield); [MW_CallLambda,MW_PutLambda]=blslambda(Price, Strike, Rate, Time, Volatility,Yield);
Note that NAG doesn’t output the elasticity (Lambda) directly but it is trivial to obtain it using values that it does output. Also note that as far as I can tell, NAG outputs more Greeks than the Financial Toolbox does.
I’m not going to show the entire output of the above program because there are a lot of numbers. However, here are the Put values as calculated by NAG shown to 4 decimal places. I have checked and they agree with the Financial Toolbox to within numerical noise.
NAG_Put =0.1350 NAG_PutDelta =-0.0419 NAG_PutLambda =-15.5066 NAG_Gamma =0.0119 NAG_Vega =2.2361 NAG_PutTheta =-1.1187 NAG_PutRho =-0.5572 NAG_PutCrho = -0.5235 NAG_PutVanna =-0.4709 NAG_PutCharm =0.2229 NAG_PutSpeed =-0.0030 NAG_PutColour =-0.0275 NAG_PutZomma =0.0688 NAG_PutVomma =20.3560