Re: Q: Quantil calculation in IDL? [message #17467 is a reply to message #17379] |
Mon, 25 October 1999 00:00   |
m218003
Messages: 56 Registered: August 1999
|
Member |
|
|
In article <7uensq$8vu$1@usenet.bham.ac.uk>,
James Tappin <sjt@star.sr.bham.ac.uk> writes:
> Joerg Mosthaf wrote:
>> Hi,
>> I have been searching the help files and David Fannings great book,
>> but I can't
>> find a way to calculate 25%- and 75%-quantils.
>> ...
Please find attached my PERCENTILES function which uses the SORT method
and allows for calculating arbitrary percentiles.
Cheers,
Martin.
[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[ [[[[[[[
[[ Dr. Martin Schultz Max-Planck-Institut fuer Meteorologie [[
[[ Bundesstr. 55, 20146 Hamburg [[
[[ phone: +49 40 41173-308 or 416 [[
[[ fax: +49 40 41173-298 [[
[[ martin.schultz@dkrz.de [[
[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[ [[[[[[[
; $Id: percentiles.pro,v 1.11 1999/05/20 16:20:42 mgs Exp $
;----------------------------------------------------------- --
;+
; NAME:
; PERCENTILES
;
; PURPOSE:
; compute percentiles of a data array
;
; CATEGORY:
; statistical function
;
; CALLING SEQUENCE:
; Y = PERCENTILES(DATA [,VALUE=value-array])
;
; INPUTS:
; DATA --> the vector containing the data
;
; KEYWORD PARAMETERS:
; VALUE --> compute specified percentiles
; default is a standard set of min, 25%, median (=50%), 75%, and max
; which can be used for box- and whisker plots.
; The values in the VALUE array must lie between 0. and 1. !
;
; OUTPUTS:
; The function returns an array with the percentile values or
; -1 if no data was passed or value contains invalid numbers.
;
; SUBROUTINES:
;
; REQUIREMENTS:
;
; NOTES:
;
; EXAMPLE:
; x = (findgen(31)-15.)*0.2 ; create sample data
; y = exp(-x^2)/3.14159 ; compute some Gauss distribution
; p = percentiles(y,value=[0.05,0.1,0.9,0.95])
; print,p
;
; IDL prints : 3.92826e-05 0.000125309 0.305829 0.318310
;
; MODIFICATION HISTORY:
; mgs, 03 Aug 1997: VERSION 1.00
; mgs, 20 Feb 1998: - improved speed and memory usage
; (after tip from Stein Vidar on newsgroup)
;
;-
; Copyright (C) 1997, Martin Schultz, Harvard University
; This software is provided as is without any warranty
; whatsoever. It may be freely used, copied or distributed
; for non-commercial purposes. This copyright notice must be
; kept with any copy of this software. If this software shall
; be used commercially or sold as part of a larger package,
; please contact the author to arrange payment.
; Bugs and comments should be directed to mgs@io.harvard.edu
; with subject "IDL routine percentiles"
;----------------------------------------------------------- --
function percentiles,data,value=value
result = -1
n = n_elements(data)
if (n le 0) then return,result ; error : data not defined
; check if speficic percentiles requested - if not: set standard
if(not keyword_set(value)) then value = [ 0., 0.25, 0.5, 0.75, 1.0 ]
; create a temporary copy of the data and sort
; tmp = data
; tmp = tmp(sort(tmp))
; NO: simply save the sorted index array
ix = sort(data)
; loop through percentile values, get indices and add to result
; This is all we need since computing percentiles is nothing more
; than counting in a sorted array.
for i=0L,n_elements(value)-1 do begin
if(value(i) lt 0. OR value(i) gt 1.) then return,-1
if(value(i) le 0.5) then ind = long(value(i)*n) $
else ind = long(value(i)*(n+1))
if (ind ge n) then ind = n-1L ; small fix for small n
; (or value eq 1.)
; if(i eq 0) then result = tmp(ind) $
; else result = [result, tmp(ind) ]
; ## change number 2
if(i eq 0) then result = data(ix(ind)) $
else result = [result, data(ix(ind)) ]
endfor
return,result
end
|
|
|