[ncl-talk] Fwd: Re: Extract upper high tail of probability distribution
Sri.durgesh Nandini-Weiss
sri.durgesh.nandini-weiss at uni-hamburg.de
Thu Dec 12 03:56:08 MST 2019
Hello everyone and sorry for the re-post,
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%), it
extracts the correct data BUT: in the attached diagram the y axis/labels
do not fit the area extracted. Would someone know how to also center
these (high tail between 95-100percentiles) onto the mean?
THe aim is to understand if its possible to see which change in
statistics is driving this high tail (variance, kurtosis or skewness?).
Below is my script for calculating and plotting this.
Would appreciate any advice.
f = addfile ("anom_hist_zo_slr.nc", "r")
hist_anom = f->hist_anom
printVarSummary(hist_anom)
f1 = addfile ("anom_rcp45_zo_slr.nc", "r")
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
be added
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)
;================================================================================
; add text
;================================================================================
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/20191212/b9190b03/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/20191212/b9190b03/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/20191212/b9190b03/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/20191212/b9190b03/attachment-0003.pdf>
-------------- next part --------------
_______________________________________________
ncl-talk mailing list
ncl-talk at ucar.edu
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
More information about the ncl-talk
mailing list