VIDEO TRANSCRIPT: Hello, this is Matt Asher from StatisticsBlog.com. I’m going to be reviewing Mathematica 9, from Wolfram Research. In particular, I’ll be focusing on using it with R and to do Monte Carlo simulations and other statistical work. You can find a full transcript of this video at my blog, including the source code and links to all of the webpages I mention.
Before I begin I’d like to thank Jeff Hara and Andy Ross from Wolfram for their time. Also thanks to the folks at the Mathematica Stack Exchange, who helped with a couple of my questions.
I’m going to get started with a blank notebook. I’m going to clear out all of the variables that exist. I’ve found sometimes that if you have existing variables that can cause problems.
ClearAll["Global`*"]
After each line I’m hitting Shift+Enter to run the command, if you just hit enter Mathematica won’t run things yet.
So I’ve cleared my variables and I’m going to run
Needs["RLink`"]
which will bring in the link between Mathematica and R.
InstallR[]
I’m going to make sure it’s installed.
REvaluate["R.Version()"]
And then I’m going to run a test command here to make sure everything is up and running. As you can see this is the version of R I’m running and the connection succeeded.
Note that the free version of Mathematica, the evaluation version, doesn’t come with support for R, so if you want to test out Mathematica’s and its interactions with R you either have to have to buy the full version or maybe if you call up or contact Wolfram they’d be willing to let you have a free evaluation version that is full and allows you to test out R.
So how does the interface between R and Mathematica work?
Basically, you can run operations in R, then store the results to variables in R. You can also pass data types back and forth between R and Mathematica.
Here I’m setting a variable and this variable is set in R, not in Mathematica
RSet["hal", 9000]
So if I were to type just
hal
There is no response back. This is Mathematica’s way of saying that the variable is undefined or that it doesn’t know what to do with your input. So to get back information from R we have to use:
REvaluate["hal"]
We are putting “hal” in quotes so we are parsing this in R and not in Mathematica.
For example we can do things like grab a dataset from R
iris = REvaluate["iris"]
I’m grabbing the famous “iris” dataset in R and I am pulling it into Mathematica.
or I could do things like evaluate a command in R:
REvaluate["det(matrix(sample(-50:49),nrow=10))"]
and bring back the results. This grabs the determinant of a random matrix.
We can even do things like create our own functions in R, and this gets put into a variable in Mathematica.
perfectSample = RFunction["function(n, dist, ...) match.fun(paste('q', dist, sep=''))((1:n) / (n+1), ...)"]
This function creates a perfect sample of the length that you specify of a particular distribution. Then we can call that function directly in Mathematica.
perfectSample[100, "pois", 10]
and the results are returned.
Of course, if we just wanted to do things in R, we would be continuing to just use R, instead of having to learn this new interface between R and Mathematica. So then what can we do in Mathematica that goes beyond what we can easily do in R?
One of the biggest advantages to using Mathematica is that you get access to tools for creating interactive plots and simulations that can be changed on the fly.
I’m going to do an example using the Benini Distribution, which, according to Wolfram’s web page, can be used to model the weight of cats.
So to do that, what I’m going to do is use the Mathematica command “Manipulate”
Manipulate[Block[{data, dist, kmd},
data = RandomVariate[dist = BeniniDistribution[\[Alpha],\[Beta], \[Sigma]], n];
kmd = KernelMixtureDistribution[data, h, MaxMixtureKernels -> All];
Plot[{PDF[dist, x], PDF[kmd, x]}, {x, 0, xRng},
PlotRange -> Automatic]], {{\[Alpha], 1}, .01, 5}, {{\[Beta], 1}, 0.01,
5}, {{\[Sigma], 1}, .01, 2}, {{n, 100}, {10, 100, 1000, 10000}}, {{xRng, 5}, 1,
10}, {{h, .5}, 0.01, 1}]
And then I get back the results and what I’ve got here is a live Monte Carlo simulation where I am specifying the different parameters of the distribution and I’m also specifying how many variates that I’m creating. This is the smoothing, the kernel bandwidth that I am adjusting.
And I can adjust the size of it here. Make it bigger. And do all of these adjustments on the fly.
As you can see, you’ve got some good power here for generating interactive plots and simulations. You can do these natively in Mathematica, or you do live manipulation of results coming from R. This example comes from the Mathematica guys:
mathematicaRPlotWrapper = RFunction["function(filename, plotfun){
pdf(filename)
plotfun()
dev.off()
}"];
Clear[getRPlot];
getRPlot[plotFun_RFunction] :=
With[{tempfile = FileNameJoin[{$TemporaryDirectory, "temp.pdf"}]},
If[FileExistsQ[tempfile], DeleteFile[tempfile]];
mathematicaRPlotWrapper[tempfile, plotFun];
If[! FileExistsQ[tempfile], Return[$Failed]];
Import[tempfile]];
myFun[t_] :=
Show[#, ImageSize -> Medium, PlotRange -> All] &@
getRPlot[RFunction["function(){
x<- seq(1," <> ToString@t <> ",0.1)
y<- sin(x)
plot(y~x)
}"]]
What's going to happen here is I am calling an R function, doing all of my calculations, bringing them back into Mathematica.
I forgot the "Manipulate" part:
Manipulate[myFun[t], {t, 2, 10}]
So here we go. And what's happening is everything is being sent to R for processing then coming all the way back to Mathematica. As you can see even though we are making that round trip the results are coming back at a good pace, it's almost instantaneous this particular example.
What about speed though, more generally?
I tried creating some random variates like I did in my examination of JavaScript versus R. So I'm going to create 10 million random variates from a Normal distribution
Timing[data = RandomVariate[NormalDistribution[], 10^7];]
and that takes about a quarter of a second, a little bit more, which is about twice as fast as R for generating those.
But then let's look at more of a worst-case scenario, a bunch of nested loops.
Timing[
l = 0;
For [i = 0, i < 10^2, i = i + 1,
For[j = 0, j < 10^2, j = j + 1,
For[k = 0, k < 10^2, k = k + 1,
l = l + 1
]
]
]
]
Three of them, a total of a million loops it's going through, and this takes about 1.2 seconds, and it will increase by an order of magnitude if I add an order of magnitude in here. That is slow, it's about twice as slow as the same code would take to run in R, a language not known for it's speed with loops. Of course, this is generally speaking not what you want to do if you want to run fast code. Mathematica itself on their website advises against using these kinds of procedural codes.
But at the same time I've found that there are times when there really is no other way to do things, especially if you are doing simulations with a bunch of objects that take place over time, that are iterative, in which case you do need a programming structure that is like this, so you'll have to keep this in mind in terms of the speed.
Beside the live graphical manipulation of the results, another key benefit to using Mathematica is that you can do direct probability calculations, and sometimes even find closed form solutions for combining random variables.
I've got some code here that determines the probability that one standard normal variable will be greater than another.
Probability[
Subscript[x, 1] > Subscript[x,
2], {Subscript[x, 1] \[Distributed] NormalDistribution[0, 1],
Subscript[x, 2] \[Distributed] NormalDistribution[0, 1]}]
And it is, of course, one-half. That's a simple example. We can do things that are more complicated. I'm going to look here for a more general solution when you have two Normal variables with means μ1 and μ2 and we're going to try and find the probability that one of them is greater than the other.
Probability[
Subscript[x, 1] > Subscript[x,
2], {Subscript[x, 1] \[Distributed]
NormalDistribution[Subscript[\[Mu], 1], 1],
Subscript[x, 2] \[Distributed]
NormalDistribution[Subscript[\[Mu], 2], 1]}]
As you can see it's taking more time to run this calculation. And eventually we get back the response. And we do have a closed form solution.
"Erfc" stands for the complementary error function.
Unfortunately, not all of the problems that I tried to do work and sometimes they just freeze up completely.
Probability[
Subscript[x, 3] > Subscript[x,
4], {Subscript[x, 3] \[Distributed] PoissonDistribution[5],
Subscript[x, 4] \[Distributed] PoissonDistribution[10]}]
Here I'm running a query to Mathematica to try and find the probability that a Poisson(5) will be greater than a Poission(10) random variable. I found that this just freezes up and will not return a response. Though I only waited a certain number of minutes. Actually, one time it locked my compter up entirely, the other time I gave up after a few minutes. I'm going to hit Alt-comma to abort the command and back out of that.
So, by comparison to R, you can do the same calculation of two Poissions. I'm going to make sure that's run in R:
x = rpois(10^6,5); y=rpois(10^6,10); z = x
(NOTE: This code actually finds the probability that a Poission(5) is LESS than a Poission(10))
As you can see I've run that a couple times already, and it takes about .9 seconds to run this. Of course, this is not trying to find a closed form solution, but for me anyway, these are perfectly good solutions, numerical solutions, to the problems.
So, going back in to Mathematica. Besides getting closed form solutions for probabilities, you can combine distributions, do convolutions and see what kind of formula you get back out.
I found this worked fairly well for simple, well-know distributions:
dist = TransformedDistribution[u + v,
{u \[Distributed] NormalDistribution[μ1, σ1],
v \[Distributed] NormalDistribution[μ2, σ2]}];
PDF[dist, x]
I do it here for the Normal. I'm adding two Normally distributed variables together and we get back out very quickly the formula for that.
But what happens if we try to work with a less well known distribution. In fact, lets go ahead and see what happens if we want to add together cats and find out the final weight distribution.
dist = TransformedDistribution[u + v,
{u \[Distributed] BeniniDistribution[\[Alpha]1,\[Beta]1,\[Sigma]1],
v \[Distributed] BeniniDistribution[\[Alpha]2, \[Beta]2,\[Sigma]2]}];
PDF[dist, x]
And, apparently, cats don't add so well. I tried this one as well a couple times and wasn't able to get results returned from Mathematica unfortunately.
So, besides these bugs and issues, there are a couple other significant downsides to Mathematica. Sharing results over the Internet can be done: you can export notebooks to HTML, but people if they want to use the interactive graphs they'll need to have a special browser plugin installed. I asked the guys from Wolfram if they know what percent of web users already have it installed. They didn't know, I suspect the number is very low, much lower than, say, the number who have Flash installed. Of course R by itself doesn't provide much support for sharing results over the web, though Rstudio makes something called Shiny that can do some exporting over the web. I haven't checked that out very much yet. I plan to do that soon.
So, beyond sharing, the biggest impediment to using Mathematica on a daily basis is the interface. The GUI. I'll bring in some of the buttons here, the palettes of buttons. Overall the look is extremely crude. Things are disordered. The floating palettes have buttons of various sizes different places, and it looks essentially like something you might have hacked together in Visual Basic 15 years ago and then hadn't been touched since. Clearly, the modern era of interface design has passed Mathematica by. At this point even open source programs that began with horrific interfaces, like the Gimp, have now begun to focus on making their interfaces look good. And looks may seem superficial, but if you are going to be working with an interface like this, if you are going to be paying $1000 for a license of Mathematica, I don't think it's too much to expect that the design be easy to use, that it be easy to find things, and that there are cues as to where to find things. Color would certainly go a long way to help with that, as would other cues to help you.
I'm going to bring back in Mathematica here.
Based on the GUI, I wonder... One more thing about the GUI, if you're moving the palettes along it doesn't do ghosting, so it pops back in.
So, despite these issues, I can see using Mathematica as an occasional tool to find exact results for probabilities and distributions, or quickly testing out how changing parameters affects a distribution's shape. For now though, I'm going to continue experimenting with JavaScript as an alternative to R that runs code quickly and also I'll be looking some more into Shiny.
Make sure to check out the links related to this video at StatisticsBlog.com, and if you like these videos click on the subscribe button.
Tags: Cat weights, Mathematica, r
I use both Mathematica and R. I’ve used Mathematica for years. Once you get used to the interface it’s not that bad. I could care less about the look of the buttons. What bothers me is that all of the creative energy from Wolfram is going into Alpha. The old-school design is meelry a symptom of the lack of attention it gets. They do the bare minimum to keep users around, like adding Rlink. They must figure everyone who uses Mathematica is locked in or they need to use it for school, so making it attractive for new users is pointless.
You can make a natural language inquiry to Mathematica and when data comes, each rectangle is called a ‘subpod’. If you right-click the subpod you can select ‘Paste as computable data’ and that is as far as Alpha gets to old-school Mathematica stuff.
Though I agree the interface needs some serious work, I hardly think that the updates to Mathematica are a bare minimum to keep users hooked. The updates to NDSolve in version 9 are massive and life changing for me. The improvements to network functionality is also huge. I’m sure different users could give there own examples. I have many issues with Mathematica, but I don’t think it is fair to say that they barely maintain it. It is clearly there flagship product that is the basis for alpha and everything else at wolfram.
I’m a Mathematica user trying to get into R because Mathematica is not the standard of my field (economics). I learned my econometric computation when economics faculties in Mexico were barely discovering SPSS and Stata, and frankly chose Mathematica because I wanted a scientific workspace that could do more than statistics.
I agree with Mark above on Alpha – although there are some ways to bring WolframAlpha code to your Mathematica notebook, I wish they would update things like CountryData instead of expecting me to get ‘disparate sources’ from Alpha to build a time series. What I do love about Mathematica is Premier Service, which for a yearly fee gives you free software upgrades and technical assistance from the guys at Wolfram.
on RLink: I am a mac user and was dismayed to find out that you can run a fully-fledged R (or your own distro) only in the Windows environment. So, is Wolfram showing us their love for R and open source through Microsoft only? Weird. I dismantled my Mathematica windows versions and frankly, I am not too enthused to do stuff that I can already do in Mathematica in R. What I would like to do is stuff that Mathematica struggles with and where R excels, such as panel data analysis, website scraping or big data.
So up until the moment where I can install packages on my RLink-ed R, RLink is pretty much an amusement but nothing more. Good for learning basic R from Mathematica but useless to actually do any work.
I should mention up front that I work for Wolfram Research. I wanted to point out a couple of things that I think were missed here.
First, I agree completely that For loops are slow in Mathematica. As you said, we generally recommend avoiding them as more elegant and fast functional solutions usually exist. When you really want procedural loops and machine precision you should consider using Compile. If you have a C compiler like Visual Studio installed Mathematica can compile directly to C.
Here is your example using compile and compiling to C. I suspect the result is faster than what you might get in R.
f = Compile[{}, Block[{l = 0, i, k, j}, For[i = 0, i < 10^2, i = i + 1, For[j = 0, j < 10^2, j = j + 1, For[k = 0, k < 10^2, k = k + 1, l = l + 1] ]]; l ] , CompilationTarget -> "C"]; In[186]:= f[] // AbsoluteTiming Out[186]= {0.008001, 1000000}You also commented on the failure of Probability for Poisson random variates. I agree that this could be improved since the underlying problem is fairly simple (in terms of special functions). Testing whether a poisson[5] > poisson[10] is equivalent to testing whether the difference is greater than zero. But this difference follows a SkellamDistribution[5,10]. In Mathematica the SurvivalFunction for a SkellamDistribution has a “closed-form” solution in terms of MarcumQ.
To be fair and compare apples-to-apples lets run a Monte-Carlo simulation to compute the same probability. This is some unnecessarily cryptic functional code but I wanted to hint at its use and save some space.
In[191]:= AbsoluteTiming[ Mean[UnitStep[#] Unitize[#]] &[#[5] - #[10]] &[ RandomVariate[PoissonDistribution[#], 10^6] &] // N] Out[191]= {0.083005, 0.074268}I think you hit the nail on the head with dynamic interactivity. Mathematica really excels at this and anything can be manipulated this way.
The issue for me ultimately isn’t why use Mathematica instead of R. I think they are both powerful and incredibly useful platforms. The question is why not use both?
For some reason it doesn’t want to pick up my CompilationTarget when I enter this code. Interested readers can check the documentation on this. http://reference.wolfram.com/mathematica/ref/CompilationTarget.html
Hi Andy and everyone,
A small nit to pick: I notice that Mathematica obtains about .074 as the probability for a Poisson(5) to be greater than a Poisson(10) random variable, both from the exact result with the Skellam distribution and from the Monte Carlo simulation; from Matt’s video, R gets about 1 – 0.880 = 0.120.
I work for Maplesoft, makers of the computer algebra system Maple, and Maple gets similar values for Monte Carlo experiments. Moreover, it can also compute the exact numbers. I tried the following in Maple 17:
with(Statistics):
X := RandomVariable(Poisson(5));
Y := RandomVariable(Poisson(10));
P := Probability(X > Y);
This returns, admittedly not exactly a closed form, but an infinite summation that is an exact expression. With a bit of manipulation (see at the end of this comment), we get this value for P:
sum(5^_k * exp(-5) * GAMMA(_k + 1, 10) / GAMMA(_k + 1)^2, _k = 0 .. infinity)
where _k is the dummy variable of the summation, the single argument GAMMA function is the gamma function you would expect, and the two argument GAMMA function is the incomplete gamma function, as explained at http://www.maplesoft.com/support/help/Maple/view.aspx?path=GAMMA . The numerical value of this summation, exact to 10 decimal places, is 0.1197937523 (as obtained by running evalf on the result).
PS: This is the manipulation I performed on P to get the nice sum:
term := op(1, P):
term := simplify(term) assuming _k :: nonnegint;
P := subsop(1 = term, P);
Wouldn’t the probability you compute be P(X >= Y) not strict inequality P(X > Y)? That explains the difference here.
Ah – you’re absolutely right! Thanks Andy – I missed that.
Erik.
Hi Andy,
I am a longtime user of Mathematica and, while it has taken time to learn, I love the program now. I have used many other data analysis packages and find Mathematica’s to be by far the most comprehensive. The problem I have is in communicating my results to others. Thus, I am trying to develop stable, reliable communications between Mathematica and R/Shiny. This question may be out of scope for this forum, but I have read through the RLink documentation, but it is incomplete. In particular, the author appears knows both packages fairly well, so he/she omits fundamental basic details like how to establish the link between two instances (Mathematica and R/Shiny). Can you point to some more complete instructions?
Thank you
Although I use Mathematica daily, it’s not for stats, so I won’t comment on that.
But I was a bit surprised that you complained about the interface. I think you were a bit quick to judge it. Usability-wise the notebook paradigm is the best I have seen for *interactive work*, and it’s just beginning to be picked up by other systems (MATLAB “cells”, ipython’s htmlnotebook, Sage). I always find the traditional command line of some other systems inconvenient and limiting.
Regarding palettes, their advantage is in customizability. The only ones I use are those that I made myself, and they usually do something more complicated than just insert symbols. I never open the built-in ones. Don’t try to use them to type special symbols, use the escape shortcuts instead, e.g. ESC dist ESC for \[Distributed]
I don’t mind the notebook paradigm as such, though for interactive work using R’s shell directly or through Rstudio works fine (and quickly, with shortcuts to repeat the last command etc). What I mind is the implementation. It’s hard to select just the right level of nesting, issues of scope are confusing and non-intuitive, and feedback is poor in terms of which clumps of code are running at any moment, which have been tweaked but not re-run.
Since you’ve been using it for a while, this may be hard. But imagine designing a similar interface from scratch, and think of all the things you would *want* to have, like breadcrumbs, tabbed notebooks, detailed running indicators, fontsize/face adjustment that isn’t circa 1994, etc.