# [ncl-talk] Extract upper high tail of probability distribution

Sri.durgesh Nandini-Weiss sri.durgesh.nandini-weiss at uni-hamburg.de
Mon Dec 9 07:05:34 MST 2019

```Hello everyone,

I wanted to show 2 plots of the distribution of my SSH anomaly data.

The first plot (Hist_RCP45_Vancouver_PDF_zoslr.pdf) shows the PDF of
Vancouver, with the 5th and 95th percentiles as vertical lines.

I have problems with the second
plot(Hist_RCP45_Vancouver_95perc_zoslr.pdf), where i just want to
extract the high end/upper tail and plot this (between 95-100%), but it
doesn't show the correct data.

Below is my script for calculating and plotting this.

hist_anom    = f->hist_anom
printVarSummary(hist_anom)

rcp45_anom   = f1->rcp45_anom
printVarSummary(rcp45_anom)

dimx = dimsizes(hist_anom)
ntim = dimx(0)           ; 240
nens = dimx(1)           ; 100
nlat = dimx(2)           ; 45
nlon = dimx(3)           ; 90
lat1=nlat
lon1=nlon
nmos = 12
nyrs   = ntim/nmos       ; 20
printVarSummary(nyrs)

;==================================================================
print ("Data read in - start calculations")

;==================================================================
; Region domain lat and lon
;==================================================================

hist_anom1=hist_anom(:,:,{48},{-125})               ;(ntim, mlon)
printVarSummary(hist_anom1)
rcp45_anom1=rcp45_anom(:,:,{48},{-125})              ;(ntim, mlon)
printVarSummary(rcp45_anom1)

;=================================================================
; Do dim_avg and dim_stddev on time and ens dimension, only take
skewness and kurtosis from stat4
;=================================================================

anom1_1d = ndtooned(hist_anom1)                     ;to carry out the
dim_stat4_n on one dimension.
printVarSummary(anom1_1d)

aveX = avg (anom1_1d)
printMinMax(aveX,0)
varstdX=stddev(anom1_1d)
printVarSummary(varstdX)

anom2_1d = ndtooned(rcp45_anom1)
printVarSummary(anom2_1d)

aveY = avg (anom2_1d)
printMinMax(aveY,0)
varstdY=stddev(anom2_1d)
printVarSummary(varstdY)

;================================================================================
; calculate 95-100 percentiles: extract mean+2standard devitations
;===============================================================================
; percentiles (requires sorting the array)

x95= aveX+2*varstdX
;no = dim_pqsort(anom1_1d,2)                      ;; sorted array is
anom1_nomiss!!
i95  =  ind(anom1_1d.ge.x95)
xHigh = anom1_1d(i95)                 ; values
printVarSummary(xHigh)

y95= aveY+2*varstdY
;it = dim_pqsort(anom2_1d,2)                      ; sorted array is
anom1_nomiss!!
i95Y  =  ind(anom2_1d.ge.y95)
YHigh =   anom2_1d(i95Y)              ; values
printVarSummary(YHigh)

;==================================================================

TTstat1 = dim_stat4(xHigh)
printVarSummary(TTstat1)

TTstat2 = dim_stat4(YHigh)
printVarSummary(TTstat2)
print("calculate PDFs=========================")

;================================================================================
; calculate empirical PDF
;================================================================================
opt          = True
opt at bin_min = -25.
opt at bin_max =  25.

pdf1 = pdfx(xHigh, 100, opt)
printVarSummary(pdf1)
pdf2 = pdfx(YHigh, 100, opt)
printVarSummary(pdf2)

;================================================================================
; prepare data for plotting
;================================================================================

nVar    = 2
nBin    = pdf1 at nbins          ; retrieve the number of bins

xx      = new ( (/nVar, nBin/), typeof(pdf1))
xx(0,:) = pdf1 at bin_center     ; assign appropriate "x" axis values
xx(1,:) = pdf2 at bin_center

yy      = new ( (/nVar, nBin/), typeof(pdf1))
yy(0,:) = (/ pdf1 /) * 0.01   ;; convert % to absolut
yy(1,:) = (/ pdf2 /) * 0.01

;================================================================================
; Plot PDFs...seperate plots
;================================================================================

plotfile="Hist_RCP45_Vancover_95perc_zoslr"
wks_type = "pdf"
wks = gsn_open_wks(wks_type,plotfile)

res                   = True
res at gsnMaximize       =True
res at vpHeightF         = 0.7                ; change aspect ratio of plot
res at vpWidthF          = 1.

res at tmYLAutoPrecision    = False          ; no auto precision
res at tmYLPrecision        = 2              ; set the precision

res at xyLineThicknesses = (/5.,5./)
res at xyDashPattern = 0               ; dashed/solid(=0)
res at trXMinF = 0             ; set minimum X-axis value
res at trXMaxF = 15             ; set maximum X-axis value
res at trYMaxF = 0.11              ; set maximum Y-axis value

res at tiXAxisFontHeightF = 0.022
res at tiYAxisFontHeightF = 0.022
res at tiYAxisString = "probability %"         ; axis string
res at tiXAxisString = "SSH anomaly [cm]"      ; axis string

res at xyLineColors      = (/"blue","red"/)
res at gsnXRefLineDashPattern =1
res at gsnXRefLine           = 0.0

res at gsnDraw           = False              ; don't draw so text can
res at gsnFrame          = False
res at gsnStringFontHeightF = 0.026
res at gsnCenterString        = "City: "+region +" (Lat:" +lat1 +",
Lon:"+lon1 +")"
plot = gsn_csm_xy (wks, xx, yy , res)

;================================================================================
;================================================================================

txres               = True
;txres at txFont        = "duplex_roman"
txres at txJust        = "CenterLeft"
txres at txFontHeightF = 0.019
;txres at txFuncCode    = "~"

;; with axis lables

txres at txFontColor = "blue"
gsn_text_ndc(wks,"1986-2006",0.19,0.645,txres)
gsn_text_ndc (wks,"Mean     = "+decimalPlaces(TTstat1(0),3,True),
0.192, 0.61,txres)
gsn_text_ndc (wks,"Variance = "+decimalPlaces(TTstat1(1),3,True),
0.192, 0.575 ,txres)
gsn_text_ndc (wks,"Skewness = "+decimalPlaces(TTstat1(2),3,True),
0.192, 0.54 ,txres)
gsn_text_ndc (wks,"Kurtosis = "+decimalPlaces(TTstat1(3),3,True),
0.192, 0.505 ,txres)

txres at txFontColor = "red"
gsn_text_ndc(wks,"RCP4.5: 1980-2100",0.698,0.645,txres)
gsn_text_ndc (wks,"Mean     = "+decimalPlaces(TTstat2(0),3,True),
0.70, 0.61,txres)
gsn_text_ndc (wks,"Variance = "+decimalPlaces(TTstat2(1),3,True),
0.70, 0.575 ,txres)
gsn_text_ndc (wks,"Skewness = "+decimalPlaces(TTstat2(2),3,True),
0.70, 0.54 ,txres)
gsn_text_ndc (wks,"Kurtosis = "+decimalPlaces(TTstat2(3),3,True),
0.70, 0.505 ,txres)

draw (plot)
frame (wks)

On 12/5/19 5:44 PM, Dennis Shea wrote:
> If the source array is 'x[*]' and the distribution in normally
> distributed with xMean and xStd, the.
>
>    x95 = xMean+2*xStd   ; scalar
>    i95  =  ind(x.ge.x95)
>    xHigh = x(i95)              ; values
>    print(xHigh)
> ===
> If 'x' is multidimensional [ 'X' ] then
>
>    x = ndtooned(X)
>
>
>
>
> On Thu, Dec 5, 2019 at 6:53 AM Sri.durgesh Nandini-Weiss via ncl-talk
> <ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>> wrote:
>
>     Hello everyone,
>
>     I am trying to extract the upper high end tail (95
>     percentiles-100percentiles) (shaded area below) of the probability
>     distribution.
>
>     Is one way: Area95-100percentiles= mean+2standard deviations?
>
>     Does anyone know how to extract (resampled) and plot this high
>     end/upper tail from a normal distribution?
>
>
>
>     --
>     Dr. Sri Nandini-Weiß
>     Research Scientist
>
>     Universität Hamburg
>     Center for Earth System Research and Sustainability (CEN)
>     Cluster of Excellence 'Climate, Climatic Change, and Society' (CLICCS)
>
>     Bundesstrasse 53, 20146 Hamburg
>     Tel: +49 (0) 40 42838 7472
>
>     _______________________________________________
>     ncl-talk mailing list
>     ncl-talk at ucar.edu <mailto:ncl-talk at ucar.edu>
>     List instructions, subscriber options, unsubscribe:
>     http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
--
Dr. Sri Nandini-Weiß
Research Scientist

Universität Hamburg
Center for Earth System Research and Sustainability (CEN)
Cluster of Excellence 'Climate, Climatic Change, and Society' (CLICCS)

Bundesstrasse 53, 20146 Hamburg
Tel: +49 (0) 40 42838 7472

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191209/73ea87f8/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: igmpcfpjcpojomll.jpeg
Type: image/jpeg
Size: 43896 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191209/73ea87f8/attachment-0001.jpeg>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Hist_RCP45_Vancouver_PDF_zoslr.pdf
Type: application/pdf
Size: 91281 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191209/73ea87f8/attachment-0002.pdf>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Hist_RCP45_Vancover_95perc_zoslr.pdf
Type: application/pdf
Size: 83519 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191209/73ea87f8/attachment-0003.pdf>
```