[Thread Prev][Thread Next][Index]

Re: [ferret_users] probability density function

Hi Izidine,
There are several scripts that do histograms in the Ferret distribution, which do this general kind of thing. (Enter "go/help frequency_histogram.jnl" for information about one of them which makes a bar chart). Searching the archives for "histogram" may bring you other ideas.

Here is another script, which is an updated version of the script called histogram.jnl. (The old histogram.jnl script does the same thing but as written it has limitations on the number of data it uses, and it writes out intermediate files which is no longer necessary with newer Ferret functions.) Save this and make a copy of it; then you can work with it and change the commands if you like to suit what you're doing. The expression can be any shape, that is it can be a multi-dimensional variable.

! histogram_pdf.jnl
\CANCEL MODE VERIFY ! based on the old histogram.jnl
! Upgraded 10/2009 for Ferret v6.3+
! Description: generate and plot a frequency histogram from a FERRET variable
! Note that internal settings are made for:
! 1) computing the PDF (see the definition of hpdf below), and
! 2) for smoothing the histogram that is plotted, which is the smoothed PDF.

! Usage:  GO histogram_pdf expression
! Example: ! USE levitus_climatology
!   GO histogram_pdf temp[k=1:2]

! check that all necessary input arguments were supplied
QUERY/IGNORE $1"<Usage:GO histogram_pdf expression"

! produce a sorted, numbered set from the data in the expression
! hsort is a sorted list of all the valid data
! hcount is a list of indices of the valid data

LET ns = `hsort[i=@ngd]`
LET hval = hsort[i=1:`ns`]

LET hcount = i[gx=hval,i=1:`ns`]

! define variables needed for the histogram
! hcdf - a normalized "counter" that increments for each data point
! hpdf - the raw probability density function for the users data
! The PDF approximates the derivative of the CDF as
! d/dx(CDF) = d/di(CDF) / d/di(HVAL) where HVAL are the sorted data values.
! An arbitrary delta i of 10 points is used.  Large data sets could use a
! larger delta i.
! e.g. LET hpdf = 50/((hval[i=@shf:+25]-hval[i=@shf:-25])*hcount[x=@max])

LET hcdf = hcount/hcount[x=@max] ! normalized for cumulative prob. dens. fct
LET hpdf = 10/((hval[i=@shf:+5]-hval[i=@shf:-5])*hcount[x=@max])

! various results can be plotted from this
! plot/vs hval,hcdf              ! cumulative probability density function
! plot/vs hval,hpdf              ! approximate probability density function
! plot/vs hval,hpdf[i=@sbx:11]   ! histogram: smoothed PDF

! For large data sets smoothing can be increased using a larger @SBX argument

LET/TITLE="Probability Density Function" vval = hpdf[i=@sbx:11]
SET VAR/TITLE="`($1),return=title` (`$1,return=units`)" hval

! Do an "underlay" plot of the original variable (but with a /VLIM setting
! so that no data points are plotted). This puts the region information in
! the upper left, and the dataset info in the upper right.

DEFINE VIEWPORT hvp1    ! just like "full"
DEFINE VIEWPORT hvp2    ! just like "full"

LET varmax = `hsort[i=@max]`

PLOT/VS/NOAX/VLIM=`varmax+1`:`varmax+2` ($1),($1)


PLOT/VS/NOLAB hval,vval

! Label the axes.
LABEL/NOUSER `($ppl$xlen)/2`,-0.8, 0, 0, 0.12, @AC`($1),return=title` (`$1,return=units`)
LABEL/NOUSER -0.8,`($ppl$ylen)/2`,0,90,0.12, @ACProbability Density Function

SAY Various results can be plotted from the definitions made by this script:

SAY PLOT/VS hval,hcdf              ! cumulative probability density function
SAY PLOT/VS hval,hpdf ! approximate probability density function
SAY PLOT/VS hval,hpdf[i=@SBX:11]   ! histogram: smoothed PDF


Izidine Pinto wrote:

Hi  Ansley

The data that I have are outputs of model simulations over southern
east africa. The 31 values are days.
What I would like to do is to produce PDF for precipitation taking in account
all domain for an interval of time (e.g 1979-1999).
I don't how to do and if its possible.

Thank you

[Thread Prev][Thread Next][Index]

Contact Us
Dept of Commerce / NOAA / OAR / PMEL / TMAP

Privacy Policy | Disclaimer | Accessibility Statement