superimposing Gaussian PDF onto histogram

From: Kerrie Geil <kgeil_at_nyahnyahspammersnyahnyah>
Date: Mon Sep 12 2011 - 18:02:44 MDT

Hi,

I'm trying to superimpose gausian and gamma PDFs on top of a histogram to
show which is a better fit for a course I'm taking. I'm plotting the
histogram and then using gsn_add_polyline to add the PDFs, but can't get the
PDFs to scale correctly over the histogram. If I plot the PDFs alone, they
look fine. Is there something with gsn_histogram that will not allow this?

My code is attached along with the small dataset I'm using.

Thanks for the help!

Kerrie

;###############################################
load "$NCARG_NCARG/nclscripts/csm/gsn_code.ncl"
load "$NCARG_NCARG/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin

;##### READ DATA ####
lees_ferry_data=asciiread("/Users/geil/Documents/stats_class/HW1/lees_ferry.txt",-1,"string")
;#### EDIT LOCATION HERE #####

nrows_lees=dimsizes(lees_ferry_data)
time=new((/nrows_lees/),string)

;##### CREATE A TIME COORDINATE #####
do i=0,nrows_lees-1
  if(toint(str_get_field(lees_ferry_data(i),2," ")).lt.10)then
    time(i)=str_get_field(lees_ferry_data(i),1,"
")+"0"+str_get_field(lees_ferry_data(i),2," ")
   else
     time(i)=str_get_field(lees_ferry_data(i),1,"
")+str_get_field(lees_ferry_data(i),2," ")
  end if
end do

;##### SEPARATE OUT THE STREAM FLOW DATA #####
flow_lees=tofloat(str_get_field(lees_ferry_data,3," "))
flow_lees!0="time"
flow_lees&time=time

;##### SEPARATE WINTER DATA #####
DJFindicesL=ind( (toint(str_get_cols(flow_lees&time,4,5)) .le.2) .or.
(toint(str_get_cols(flow_lees&time,4,5)).eq.12) )
DJF_L=flow_lees(DJFindicesL)

;##### CALCULATE GAUSSIAN #####
mean_DJFL=sum(DJF_L)/dimsizes(.not.(ismissing(DJF_L)))
stddev_DJFL=(sum((DJF_L-mean_DJFL)^2)/num(.not.(ismissing(DJF_L))))^0.5

pi=3.14159
qsort(DJF_L)

x1=ispan(0,30000,100)
gausDJFL=(1./(stddev_DJFL*((2*pi)^0.5)))*exp(-(((x1-mean_DJFL)^2)/(2*stddev_DJFL^2)))

;#################################################################################
;##### PLOTTING SECTION
##########################################################
;#################################################################################

res=True
res@gsnMaximize=True
res@gsnHistogramSelectNiceIntervals=False
res@tmXBLabelAngleF=270.
res@gsFillColor=(/"transparent"/)
res@gsEdgeColor="black"
res@tmXBLabelFontHeightF=0.012
res@tmXTBorderOn=False
res@tmYRBorderOn=False
;res@tmYRMinorOn=False
res@tmYROn=False
res@tmXTOn=False
res@tmXBLabelJust="CenterCenter"
res@tmXBMajorLengthF=0.01
res@tmYLMinorLengthF=0.01
res@tiMainFontHeightF=0.012
res@tiXAxisFontHeightF=0.012
res@tiYAxisFontHeightF=0.012

wks=gsn_open_wks("x11","DJF_L")

res@gsnDraw=False
res@gsnFrame=False

txres=True
txres@txJust="BottomCenter"
txres@txFontHeightF=0.013

plotmin=min(DJF_L)
plotmax=max(DJF_L)
h=1633.5 ;##### BIN WIDTH #####
res@gsnHistogramNumberOfBins=round((plotmax-plotmin)/h,3)
plot1=gsn_histogram(wks,decimalPlaces(DJF_L,1,True),res)
labels=""+plot1@NumInBins
dum1=gsn_add_text(wks,plot1,labels,plot1@MidBarLocs,plot1@NumInBins
+0.5,txres)

gres=True
gres@gsLineColor="blue"
mkres=True
mkres@gsMarkerIndex=17
mkres@gsMarkerSizeF=0.035
mkres@gsMarkerColor="blue"

x=fspan(0,30000,301)
dum2=new(dimsizes(gausDJFL)-1,graphic)
dum3=new(dimsizes(gausDJFL)-1,graphic)
gausDJFL=gausDJFL*102910.5 ;##### 102910.5 is No.Data Values * Bin
Width=scaling factor #####

do i=0,dimsizes(gausDJFL)-2
xvalues=(/x(i),x(i+1)/)
yvalues=(/gausDJFL(i),gausDJFL(i+1)/)
dum2(i)=gsn_add_polyline(wks,plot1,xvalues,yvalues,gres)
dum3(i)=gsn_add_polymarker(wks,plot1,x,gausDJFL,mkres)
end do

draw(plot1)
frame(wks)

plot=gsn_xy(wks,x,gausDJFL,res)
draw(plot)
frame(wks)

end
;##########################################################

-- 
Kerrie Geil
Master's Student
Department of Atmospheric Sciences
University of Arizona
PAS Building Rm 526
1118 E 4th Street
Tucson, AZ 85721



_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk

Received on Mon Sep 12 18:02:53 2011

This archive was generated by hypermail 2.1.8 : Fri Sep 16 2011 - 11:24:25 MDT