Plotting tens of thousands of functions together on the same axes

Question

I'm trying to plot a large number of functions together on the same axes. Each of these functions is the result from one iteration of a Monte Carlo simulation. Essentially, I'm varying parameters of a differential algebraic equation system and simulating the sensitivity of the solution to its parameters. The relevant piece of Mathematica code for the plot is:

Plot[
  Evaluate[{CaPlasma[t], CvPlasma[t]} /. monteCarloSolutions], 
  {t, 0, 300}, ImageSize -> Large, 
  PlotRange -> All, PlotStyle -> twoPlotColors
] // AbsoluteTiming

monteCarloSolutions in the above code is in the form of a 10,000 by 44 list of replacement rules like:

{{CaPlasma->InterpolatingFunction[...], CvPlasma->InterpolatingFunction[...], ...},
 ...    
 {CaPlasma->InterpolatingFunction[...], CvPlasma->InterpolatingFunction[...], ...}}

Below is some example output of a subset of the results (3000 iterations of the Monte Carlo simulation). I have not yet been successful in plotting all 10,000 of the iterations without crashing the Mathematica kernel (at least I think that is what happens; I don't actually get an error message). Plotting 1000 iterations took 5 minutes, 2000 took 20 minutes, and 3000 took 45 minutes.

So, basically, I'm running up against a performance issue trying to plot this many functions together in one graph. However, I also suspect that maybe this is a symptom of greater issue. Perhaps the data has outgrown this particular type of visualization and it would be better to take a different approach. I have previously tried computing summary statistics over each group of functions and then plotting a line for the mean(t) and shaded areas for ± stddev(t) (or quantiles, etc.). However, that assumes that the behavior of a group of functions is unimodal, which in some cases does not appear to be a valid assumption. I think there's probably a better way to visualize this data, but I am at a loss for identifying what type of plot would be appropriate.

Ultimately, I guess I really have a couple of questions here.

Is there a way around the performance bottleneck to make Mathematica plot tens of thousands of functions together on the same axes?
Is that even the right thing to do? If not, what alternative type of graph would better serve the purpose of visualizing this kind of data?

Evaluate the functions only at multiples of 25, and produce a DistributionChart for each point. — Alan, Dec 17 '19 at 19:56
It's not clear to me how DistributionChart could be used for my data since each Monte Carlo result is a time series and there doesn't appear to be a time axis on the DistributionChart plots. — Matt, Dec 17 '19 at 20:07
There is no point in trying to plot tens of thousands of functions; no reasonable reproduction system will have the resolution to actually show them anyway. You definitely need to find an alternative type of visualization. What features are you looking to emphasize? Would a max-min value for each function be an interesting output? In any case, you need to find a statistic to summarize your data; I'm afraid that presenting the raw runs would be meaningless for the purposes of visualization. — MarcoB, Dec 17 '19 at 20:50
@Matt The time (only each 25th periods) will be the category for the Distribution Charts. — Alan, Dec 17 '19 at 21:00
The reality is, there's a threshold where plotting past a certain number is meaningless. If you have 1,000,000 equally distributed points, plotting 10,000 of them probably gets the point across for the sake of the visualization. If you have 10,000 equally distributed lines, plotting 1,000 of them probably gets the point across. That would be my main recommendation in situations like this (that I follow all the time). Note that "equally distributed" is essential for this, though. — ktm, Dec 17 '19 at 21:00
Not sure what you mean by "However, that assumes that the behavior of a group of functions is unimodal, which in some cases does not appear to be a valid assumption.". Certainly finding a variety of percentiles for each value on the horizontal axis does not imply a unimodal distribution. In fact, you want to use percentiles when a distribution is multimodal. Depending on what feature you want to highlight (which you haven't stated) you might want to fit a SmoothKernelDistribution to each value along the horizontal axis. Then a 3D plot of all of those distributions. — JimB, Dec 17 '19 at 23:10
Separate each run in $y$, to get a large number of independent plots. (Avoid stacking them, as you have now.) . Also, use a different color for each line. — David G. Stork, Dec 18 '19 at 01:37

Mr.Wizard · Answer 1 · 2019-12-18T14:03:17.877

It's probably best to build this up in raster. Borrowing JimB's sample data:

SeedRandom[1]
z = RandomVariate[UniformDistribution[{2, 10}], 1000];
pdf = (PDF[GammaDistribution[2, #], x][[1, 1, 1]] &) /@ z;

I make a plotting function (fn) with a fixed PlotRange and reduced Opacity, then individually* Rasterize and Overlay with Fold.

*Update: plotting blocks of lines is much faster than each one alone.

fn = Plot[#, {x, 0, 20}, AxesLabel -> {"x", "y"}, PlotRange -> {0, 0.2}, 
       PlotStyle -> Opacity[0.1, ColorData[97, 1]], ImageSize -> 500] &;

(* 25 lines at a time *)
Fold[Rasterize @ Overlay[{#, fn[#2]}] &, fn[-1], Partition[pdf, 25]]

Notes:

This render took under 8 seconds on my circa 2011 PC in Mathematica 10.1.
fn[-1] just gives me a base for the first overlay; a horizontal line at -1 is out of the plot range and invisible.
The aliasing seen on the upper right surface of the envelope could be eliminated by rasterizing at a larger size, then downsampling the final result. This will also smooth the font. (Use a larger font size as needed.)
Alternatively we could render the axes only once, using Axes -> False for the rest, but doing it this way was simpler and avoided alignment errors.

Is it possible to draw really thin lines in Mathematica?

score 3 · Accepted Answer · answered Dec 18 '19 at 01:30

If the functions are selected according to some random selection of the parameters, then one could look at the distribution of the function values for various values of the horizontal variable. (This might not be what you want but your request current just asks for a way to visualize the data. More specificity would be helpful.)

(* Generate a random set of 1,000 curves *)
z = RandomVariate[UniformDistribution[{2, 10}], 1000];
pdf = (PDF[GammaDistribution[2, #], x][[1, 1, 1]] &) /@ z;

(* Plot the first 20 curves *)
Plot[pdf[[1 ;; 20]], {x, 0, 20}, AxesLabel -> {"x", "y"}]

(* Find the distribution of the values of the function (y) for each value of x *)
skd = SmoothKernelDistribution[pdf /. x -> #, Automatic,
  {"Bounded", {0, ∞}, "Gaussian"}] & /@ (Range[1, 100]/5);
data = Flatten[Table[{x/5, y, PDF[skd[[x]], y]}, {x, 1, 100}, {y, 0, 0.2, 0.0005}], 1];
ListPlot3D[data, PlotRange -> All, AxesLabel -> {"x", "y", "Density"}]

Plotting tens of thousands of functions together on the same axes

2 Answers2