Moments from an Empirical Distribution

Function calc_DistributionMoment(wave wdistr, variable n [variable xmax, wave wdistrX])
 
    variable rtnv, scaled=1, limited=0, pmax
    
    if (!ParamIsDefault(wdistrX))
        scaled = 0
    endif
    if (!ParamIsDefault(xmax))
        limited=1
        if (scaled)
            pmax = x2pnt(wdistr,xmax)
        else
            pmax = x2pnt(wdistrX,xmax)
        endif
    else
        pmax = numpnts(wdistr) - 1
        if (scaled)
            xmax = pnt2x(wdistr,pmax)
        endif
    endif
    
    switch(n)
        case -2:    // rms
            duplicate/FREE wdistr, iwave, iMwave
            if (scaled)
                integrate/METH=1 iwave
                iMwave = wdistr*x^2
                integrate/METH=1 iMwave
            else
                integrate/METH=1 iwave /X=wdistrX
                iMwave = wdistr*wdistrX^2
                integrate/METH=1 iMwave /X=wdistrX
            endif
            rtnv = sqrt(iMwave[pmax]/iwave[pmax])
            break       
        case 0: // most probable
            if (scaled)
                if (limited)
                    WaveStats/M=1/Q/R=(0,xmax) wdistr
                else
                    WaveStats/M=2/Q wdistr
                endif
                rtnv = v_maxloc
            else
                if (limited)
                    WaveStats/P/M=2/Q/R=[0,pmax] wdistr
                else
                    WaveStats/P/M=2/Q wdistr
                endif
                rtnv = wdistrX[v_maxRowLoc]
            endif
            break
        default:    // all others (1-mean, 2-second, 3...)
            duplicate/FREE wdistr, iwave, iMwave
            if (scaled)
                iwave = wdistr*x^(n-1)
                integrate/METH=1 iwave
                iMwave = wdistr*x^n
                integrate/METH=1 iMwave
            else
                iwave = wdistr*wdistrX^(n-1)
                integrate/METH=1 iwave /X=wdistrX
                iMwave = wdistr*wdistrX^n
                integrate/METH=1 iMwave /X=wdistrX
            endif
            rtnv = iMwave[pmax]/iwave[pmax]
            break
    endswitch
    
    return rtnv
end