Calculating the PDF of a waveform from its samples
What I'd go with is essentially Jason R's "random resampler", which in turn is a presampled-signal based implementation of yoda's stochastic sampling.

I've used simple cubic interpolation to one random point between each two samples. For a primitive synth sound (decaying from a saturated non-bandlimited square-like signal +even harmonics to a sine) it looks like this:



Let's compare it to a higher-sampled version,



and the weird one with the same samplerate but no interpolation.



Notable artifact of this method is the overshoot in the square-like domain, but this is actually what the PDF of the sinc-filtered signal (as I said, my signal is not bandlimited) would also look like and represents the perceived loudness much better than the peaks, if this were an audio signal.

Code (Haskell):

cubInterpolate vll vl v vr vrr vrrr x
    = v*lSpline x + vr*rSpline x
      + ((vr-vl) - (vrr-vll)/4)*ldSpline x
      + ((vrr-v) - (vrrr-vl)/4)*rdSpline x
     where lSpline x = rSpline (1-x)
           rSpline x = x*x * (3-2*x)
           ldSpline x = x * (1 + x*(x-2))
           rdSpline x = -ldSpline (1-x)

                   --  rand list   IN samples  OUT samples
stochasticAntiAlias :: [Double] -&gt; [Double] -&gt; [Double]
stochasticAntiAlias rs (lsll:lsl:lsc:lsr:lsrr:[]) = []
stochasticAntiAlias (r:rLst) (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)
    = ( cubInterpolate lsll lsl lsc lsr lsrr lsrrr r )
          : stochasticAntiAlias rLst (lsll:lsl:lsc:lsr:lsrr:lsrrr:t)


rand list is a list of random variables in range [0,1].