# Re: computing bivariate frequencies

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Tue, 8 Aug 2006 09:49:16 -0600 (MDT)

> I'm trying to compute the joint pdf of two variables C and V stored in arrays
> sized (ntime x nlat x nlon). In order to obtain the absolute frequencies, I
> need to count the number of grid cells that fall within arbitralily defined
> bins for each variable. Right now I'm attempting this using a nested do loop,
> like:
>
> do ic=0,nbin_c-1
> do iv=0,nbin_v-1
> count_c_v(ic,iv)=num(ind(c.ge.c_bin(ic).and.c.lt.c_bin(ic+1).and.
> v.ge.v_bin(iv).and.v.lt.v_bin(iv+1)))
> end do
> end do
>
> I was wondering if anybody would know how to code this more efficiently on
> ncl, in order to avoid the nested loop.
>
> ----------------------------------------------------------------

[1] I do not know why you have the "ind" function
within in the "num" function. Not needed. Extra
operations so slower.

[2] I can not think of an alternative to the double do loop.
However, I did run the following: 10000 points and 2000 bins
and it ran pretty quickly ... even with 'if' and 'print'
within the inner most loop.

I deliberately chose c_bin and v_bin such that there would
be outliers. Thus, the value of 'K' will not add up to 10000.

begin
N = 10000
c = random_normal ( 0,5,N)
v = random_normal ( 0,5,N)
;print("c="+c+" v="+v)

nc_bin = 51 ; => 50 bins
nv_bin = 41 ; => 40 bins
c_bin = fspan(-10,10,nc_bin)
v_bin = fspan(-5 ,20,nv_bin)

count = new ( (/nc_bin,nv_bin/) , integer)
count = 0 ; initialize to 0

do ic=0,nc_bin-2
do iv=0,nv_bin-2
count(ic,iv) = num(c.ge.c_bin(ic) .and. c.lt.c_bin(ic+1) .and. \
v.ge.v_bin(iv) .and. v.lt.v_bin(iv+1) )
if (count(ic,iv).gt.0) then
print("ic="+ic+" iv="+iv+" count="+count(ic,iv) )
end if
end do
end do

fpc = (100.*count)/N ; frequency [%] (nc_bin,nv_bin)
printVarSummary(fpc)

K = sum( count ) ; < N
F = sum( fpc )
print("------")
print("K="+K+" F[%]="+F)
print("------")

; extra counting stuff

cmin = min(c)
cmax = max(c)
cmin = min(v)
cmax = max(v)

cbmin = min(c_bin)
cbmax = max(c_bin)
vbmin = min(v_bin)
vbmax = max(v_bin)
; count corner cases
k1 = num (c.lt.cbmin .and. v.lt.vbmin)
k2 = num (c.lt.cbmin .and. v.ge.vbmax)
k3 = num (c.ge.cbmax .and. v.lt.vbmin)
k4 = num (c.ge.cbmax .and. v.ge.vbmax)
; one outlier
k5 = num (c.lt.cbmin .and. v.ge.vbmin .and. v.lt.vbmax)
k6 = num (c.ge.cbmax .and. v.ge.vbmin .and. v.lt.vbmax)
k7 = num (v.lt.vbmin .and. c.ge.cbmin .and. c.lt.cbmax)
k8 = num (v.ge.vbmax .and. c.ge.cbmin .and. c.lt.cbmax)

KTOT = K + k1+ k2+ k3+ k4+ k5+ k6+ k7+ k8 ; should = N

print("k1="+k1)
print("k2="+k2)
print("k3="+k3)
print("k4="+k4)
print("k5="+k5)
print("k6="+k6)
print("k7="+k7)
print("k8="+k8)

print("------")
print("KTOT="+KTOT)
print("------")

if (KTOT.eq.N) then
print (" Lucky!!! All numbers accounted for!! :-)")
end if
end
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Aug 08 2006 - 09:49:16 MDT

This archive was generated by hypermail 2.2.0 : Wed Aug 09 2006 - 08:25:35 MDT