Violin plots

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.
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
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.
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 violin2.jpg
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
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.