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.

1. Is there a way around the performance bottleneck to make Mathematica plot tens of thousands of functions together on the same axes?
2. 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?

Answer 1

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"}]

Answer 2

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.

Related:

• Is it possible to draw really thin lines in Mathematica?