Violin plots

Share icon

Is it possible to make a violin plot in IgorPro? These plots show the kernel density estimation for a data distribution and can be superimposed on a box plot. The advantage to these plots is that they reveal multimodal distributions that are masked by a box plot.
An example is here in matplotlib
http://pyinsci.blogspot.co.uk/2009/09/violin-plot-with-matplotlib.html
Unless I'm missing something, I can't see an implementation in IgorPro.

My initial thought to do this is: make a histogram from each 1D wave, smooth, make a copy, multiply one version by 0.5 and the other by -0.5. There are several problems with this: dealing with multiple waves will be difficult, it will be in the wrong orientation, smoothing needs to be done relative to bin size (I think). I wondered if anyone has tackled this already? Perhaps as part of implementing the Scatter Dot Plot graphing package?

The original reference for Violin Plots is Hintze & Nelson (1998) The American Statistician 52(2):181-4.
Between DrawPoly and DrawBezier, you can pretty much draw anything you can think of.

Be sure to use axis coordinates with SetDrawEnv.


--Jim Prouty
Software Engineer, WaveMetrics, Inc.

November 24, 2014 at 08:32 am - Permalink

If you can come up with suitable waves for the halves of the violin, perhaps plot each one as a trace (you'll have to create a suitable constant Y wave to set the position of the traces). When you add a violin trace to the graph, use AppendToGraph/VERT to make it stand on its head.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

November 26, 2014 at 09:20 am - Permalink

Thank you John and Jim for the replies. I've explored the drawpoly idea, which I like because adding colour is straightforward. Once I figure out how the smoothing is done in other programs I should be able to do the violin plot using either method.

Edit: in case anyone is following this thread. I have found a very nice kernel density estimation on IgorExchange written by bech http://www.igorexchange.com/node/4749 I'll try and work up a solution to this and will post if I get something that will be useful to others.

November 30, 2014 at 11:43 pm - Permalink

I'm close to getting the violin plots working in Igor, thanks to bech's code for 1D KDE. However, I have a problem that I can't resolve and was hoping someone can help - I didn't know whether to start a new topic.

Attached are two images showing my violin plots overlaid on violin plots from R. The data are oldfaithful$waiting and oldfaithful$eruptions. My plots are in orange/red and the R plots are in black/blue. I've scaled them so you can see the difference between them. I know that I've not trimmed my violins yet.

My problem is the smoothing parameter (the standard deviation of the kernels that are summed to give the KDE). I calculate this using Silverman's rule of thumb: h= (4/3n)^1/5 * sigma. This is what R is using as the optimum smoothing parameter. Plugging this into my code (variable bw, below), I have to use h*sqrt(2) in the /WDTH flag of FastGaussTransform. This gives an "undersmoothed" violin. I can't figure out where the discrepancy comes from. FastGaussTransform definitely uses /WDTH correctly and I can't find a problem with the /RX flag (although bech's code notes that this may need attention). For the oldfaithful$waiting data if I use bw=10 in my code, Igor produces a violin that looks like the "real thing", whereas the optimum is 6.6, i.e. 4.7*sqrt(2). Any help would be much appreciated!

Function kde(w,bw,xpos)                         // 1d kernel density estimation (Gaussian kernel)
    wave w
    variable bw
    variable xpos
    
    variable n=numpnts(w), kde_norm,  xmin=wavemin(w)-2*bw, xmax=wavemax(w)+2*bw
    variable nx = min(round(100*(xmax-xmin)/bw),100)
    make /d/free/n=(n) wweights=1
    FastGaussTransform /WDTH=(bw)/RX=(4*bw) /OUT1={xmin,nx,xmax} w,wweights     // you may need to tweak /RX flag value
    wave M_FGT;  kde_norm=sum(M_FGT)*deltax(M_FGT);     M_FGT /= kde_norm
    string wnL=NameofWave(w)+"_kdeL";   duplicate /d/o M_FGT $wnL
    string wnR=NameofWave(w)+"_kdeR";   duplicate /d/o M_FGT $wnR
    string wnLx=NameofWave(w)+"_kdeLx"; duplicate /d/o M_FGT $wnLx
    string wnRx=NameofWave(w)+"_kdeRx"; duplicate /d/o M_FGT $wnRx
    string vwny=NameofWave(w)+"_kdevy"
    string vwnx=NameofWave(w)+"_kdevx"
    killwaves M_FGT
    wave wnL1=$wnL
    wave wnR1=$wnR
    wave wnLx1=$wnLx
    wave wnRx1=$wnRx
    wnLx1 =x
    wnRx1 =x
    wnR1 *=-1
    Reverse wnR1,wnRx1
    Concatenate /O /NP {wnL1,wnR1}, $vwny
    Concatenate /O /NP {wnLx1,wnRx1}, $vwnx
    wave vwny1=$vwny
    wave vwnx1=$vwnx
    setdrawenv ycoord=left, xcoord=bottom, fillpat=0
    DrawPoly vwny1[0],xpos+vwnx1[0],1,1, vwny1,vwnx1
End
violin1.jpg (43.65 KB) violin2.jpg (42.19 KB)

December 5, 2014 at 12:15 am - Permalink

I think I've solved it. The /TET flag needs including and setting to something big, e.g. 10 to increase the number of terms for Taylor expansion of Gauss. With no flag, small n examples were plotted fine, but larger waves became problematic. I'll try and get something generally useful onto Code Snippets soon.

UPDATE: code for violin plots is on code snippets http://www.igorexchange.com/node/6181

December 16, 2014 at 03:07 am - Permalink

I'm glad to see you use the FastGaussTransform. FWIW, you may want to try StatsKDE when you get your hands on IP7.

A.G.
WaveMetrics, Inc.

December 16, 2014 at 01:42 pm - Permalink