#################   this R code shows how a Chisqr distrbution 
##################   really equal to the sum of sqared normal distribution
####   see Davis for an explanation
help(Chisquare)
#####################################################
#####################################################
#####################################################
M = 10000
zip = rep(0, M)

k = 2
for(i in 1:length(zip))
{
zee = rnorm(k)
zip[i] = sum(zee^2)
}
RK = rchisq(M, k)

DZ = density(zip)
DK = density(RK)

xlim = range(c(RK, zip))
dev.set(2)
par(mfrow=c(2,1))
hist(zip, breaks=20, xlim=xlim, freq=FALSE)
lines(DZ$x, DZ$y, , col='red')
hist(RK, breaks=20, xlim=xlim, freq=FALSE)
lines(DK$x, DK$y, col='blue')


dev.set(3)
par(mfrow=c(2,1))
qqplot(zip, RK)

plot(c(DZ$x, DK$x), c(DZ$y, DK$y), type='n')
lines(DZ$x, DZ$y, , col='red')
lines(DK$x, DK$y, col='blue')


##########################################
##########################################
##########################################

xkm = runif(100, 0, 25)
ykm = runif(100, 0, 25)


XYLOC = cbind(xkm, ykm)

write.table( XYLOC , file="D:/LEES/CLASSES/GEOL_520 Data_Analysis/LAB_NOTES/xyDATA1.txt", row.names = FALSE, col.names = FALSE, quote=FALSE)



xdiv = 5
ydiv =5

plot(xkm, ykm, type='n', xlim=c(0,25), ylim=c(0,25))
abline(v=seq(from=0, to=25+1, by=xdiv), h=seq(from=0, to=25+1, by=ydiv), lty=2)
points(xkm, ykm)

N = length(xkm)
k = 0

KOUNT = rep(0, length=xdiv*ydiv)

for(i in 1:xdiv)
{
ex = (i-1)*xdiv+.5*xdiv
for(j in 1:ydiv)
{
why = (j-1)*ydiv+.5*ydiv
points(ex, why, col="red", pch=6)
k = k + 1
text(ex, why, k, pos=3)


KOUNT[k] = length(which(xkm>(ex-.5*xdiv)& xkm<=(ex+.5*xdiv) & ykm>(why-0.5*ydiv) & ykm<=(why+.5*ydiv) ))
}

}



pred = rep(N/(xdiv*ydiv),  length=xdiv*ydiv)

chi = sum(  (KOUNT-pred)^2 / pred )

dof = length(pred)

qchisq(.95, dof)


###################################


xkm = rnorm(1000, 12, 8 )
ykm = rnorm(1000, 12, 8)

points(xkm, ykm)

x = xkm[xkm>0&xkm<25]
xkm = x[1:100]

y = ykm[ykm>0&ykm<25]
ykm = y[1:100]

XYLOC = cbind(xkm, ykm)

write.table( XYLOC , file="D:/LEES/CLASSES/GEOL_520 Data_Analysis/LAB_NOTES/xyDATA2.txt", row.names = FALSE, col.names = FALSE, quote=FALSE)



##################


T1 = scan(file="D:/LEES/CLASSES/GEOL_520 Data_Analysis/LAB_NOTES/xyDATA2.txt", list(x=0, y=0))

xkm = T1$x
ykm = T1$y


plot(xkm, ykm, type='n', xlim=c(0,25), ylim=c(0,25))
abline(v=seq(from=0, to=25+1, by=xdiv), h=seq(from=0, to=25+1, by=ydiv), lty=2)
points(xkm, ykm)

N = length(xkm)
k = 0

KOUNT = rep(0, length=xdiv*ydiv)

for(i in 1:xdiv)
{
ex = (i-1)*xdiv+.5*xdiv
for(j in 1:ydiv)
{
why = (j-1)*ydiv+.5*ydiv
points(ex, why, col="red", pch=6)
k = k + 1
text(ex, why, k, pos=3)


KOUNT[k] = length(which(xkm>(ex-.5*xdiv)& xkm<=(ex+.5*xdiv) & ykm>(why-0.5*ydiv) & ykm<=(why+.5*ydiv) ))
}

}



pred = rep(N/(xdiv*ydiv),  length=xdiv*ydiv)

chi = sum(  (KOUNT-pred)^2 / pred )

dof = length(pred)

tchi = qchisq(.95, dof)

print(paste(sep = " ", chi, "compare to TEST CHI=", tchi))