library(emdbook)
library(bbmle)

setwd("...")

powderDat <- dget("powder_data.txt")
names(powderDat)

orange1 <- powderDat[[1]]
orange2 <- powderDat[[2]]
pink1 <- powderDat[[3]]
pink2 <- powderDat[[4]]
yellow1 <- powderDat[[5]]
yellow2 <- powderDat[[6]]

coords <- read.csv(".../plotmap.csv")
trees <- coords[c(2,5,6),]

orangeSourceXY <- trees[1,2:3]
yellowSourceXY <- trees[2,2:3]
pinkSourceXY <- trees[3,2:3]

orange1 <- cbind(orange1[,c(3:4,2)], sqrt((orangeSourceXY$X-orange1$X)^2 +
                                          (orangeSourceXY$Y-orange1$Y)^2))
names(orange1) <- c("X", "Y", "area", "dists")

yellow1 <- cbind(yellow1[,c(2:3,1)], sqrt((yellowSourceXY$X-yellow1$X)^2 +
                                          (yellowSourceXY$Y-yellow1$Y)^2))
names(yellow1) <- c("X", "Y", "area", "dists")

pink1 <- cbind(pink1[,c(2:3,1)], sqrt((pinkSourceXY$X-pink1$X)^2 +
                                      (pinkSourceXY$Y-pink1$Y)^2))
names(pink1) <- c("X", "Y", "area", "dists")

orange2 <- cbind(orange2[,c(3:4,2)], sqrt((orangeSourceXY$X-orange2$X)^2 +
                                          (orangeSourceXY$Y-orange2$Y)^2))
names(orange2) <- c("X", "Y", "area", "dists")

yellow2 <- cbind(yellow2[,c(3:4,2)], sqrt((yellowSourceXY$X-yellow2$X)^2 +
                                          (yellowSourceXY$Y-yellow2$Y)^2))
names(yellow2) <- c("X", "Y", "area", "dists")

pink2 <- cbind(pink2[,c(3:4,2)], sqrt((pinkSourceXY$X-pink2$X)^2 +
                                      (pinkSourceXY$Y-pink2$Y)^2))
names(pink2) <- c("X", "Y", "area", "dists")

maxdist <- max(orange1$dists, yellow1$dists, pink1$dists, orange2$dists,
               yellow2$dists, pink2$dists)

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

#negative exponential
NE <- function(x,a){
  (1/a)*exp(-x/a)}

NLL.NE <- function(x,a){
  kx <- NE(dists,a)
  -sum(log(kx))}
  
#Weibull
NLL.weib <- function(x, a, b){
  -sum(dweibull(dists, a, b, log=TRUE))
}

#log-normal
NLL.lnorm <- function(x,a,b){
  -sum(dlnorm(dists, a, b, log=TRUE))
}

#1D student-t
st <- function(x,a,b){
  (2/a)*(gamma((b+1)/2)/((b*pi)^0.5*gamma(b/2)))*(1+(x/a)^(2/b))^(-(b+1)/2)}

NLL.st <- function(x,a,b){
 kx <- st(dists,a,b)
  -sum(log(kx))}

#WALD
wald <- function(x,a,b){
  ((a/(2*pi*x^3))^0.5)*exp(-(a*(x-b)^2)/(2*b^2*x))}

NLL.wald <- function(a,b){
  kx <- wald(dists,a,b)
  -sum(log(kx))}

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

pinkMod <- NULL
orangeMod <- NULL
yellowMod <- NULL

pink <- pink1[,3:4] 
orange <- orange1[,3:4] 
yellow <- yellow1[,3:4]

##########################
for (i in 1:nrow(pink)){
  dists <- pink[i,2]
  reps <- round(pink[i,1]/10)
  pinkVec <- rep(dists, reps)
  pinkMod <- c(pinkMod, pinkVec)
}

dists <- pinkMod

mle.NE <- mle2(NLL.NE, start=list(a=10), method="BFGS")
ne.disp <- coef(mle.NE)["a"]

stud.t <- mle2(NLL.st, start=list(a=10,b=1), method="BFGS")
st.scale <- coef(stud.t)["a"]
st.shape <- coef(stud.t)["b"]
  
mle.weib <-  mle2(NLL.weib, start=list(a=10, b=1), method="BFGS")
weib.a <- coef(mle.weib)["a"]
weib.b <- coef(mle.weib)["b"]

mle.lnorm <-  mle2(NLL.lnorm, start=list(a=2, b=1), method="BFGS")
lnorm.a <- coef(mle.lnorm)["a"]
lnorm.b <- coef(mle.lnorm)["b"]

mle.wald <- mle2(NLL.wald, start=list(a=10, b=10), method="BFGS")
wald.a <- coef(mle.wald)["a"]
wald.b <- coef(mle.wald)["b"]

print(AICtab(mle.NE, mle.weib, mle.wald, mle.lnorm, stud.t,
delta=TRUE, weights=TRUE, sort=TRUE))

mle.pink1 <- mle.lnorm

##########################
for (i in 1:nrow(yellow)){
  dists <- yellow[i,2]
  reps <- round(yellow[i,1]/10)
  yellowVec <- rep(dists, reps)
  yellowMod <- c(yellowMod, yellowVec)
}

dists <- yellowMod

mle.NE <- mle2(NLL.NE, start=list(a=10), method="BFGS")
ne.disp <- coef(mle.NE)["a"]

stud.t <- mle2(NLL.st, start=list(a=10,b=1), method="BFGS")
st.scale <- coef(stud.t)["a"]
st.shape <- coef(stud.t)["b"]
  
mle.weib <-  mle2(NLL.weib, start=list(a=10, b=1), method="BFGS")
weib.a <- coef(mle.weib)["a"]
weib.b <- coef(mle.weib)["b"]

mle.lnorm <-  mle2(NLL.lnorm, start=list(a=2, b=0.5), method="BFGS")
lnorm.a <- coef(mle.lnorm)["a"]
lnorm.b <- coef(mle.lnorm)["b"]

mle.wald <- mle2(NLL.wald, start=list(a=10, b=10), method="BFGS")
wald.a <- coef(mle.wald)["a"]
wald.b <- coef(mle.wald)["b"]

print(AICtab(mle.NE, mle.weib, mle.wald, mle.lnorm, stud.t,
delta=TRUE, weights=TRUE, sort=TRUE))

mle.yellow1 <- mle.lnorm

##########################
for (i in 1:nrow(orange)){
  dists <- orange[i,2]
  reps <- round(orange[i,1])
  orangeVec <- rep(dists, reps)
  orangeMod <- c(orangeMod, orangeVec)
}

dists <- orangeMod

mle.NE <- mle2(NLL.NE, start=list(a=10), method="BFGS")
ne.disp <- coef(mle.NE)["a"]

stud.t <- mle2(NLL.st, start=list(a=10,b=1), method="BFGS")
st.scale <- coef(stud.t)["a"]
st.shape <- coef(stud.t)["b"]
  
mle.weib <-  mle2(NLL.weib, start=list(a=10, b=1), method="BFGS")
weib.a <- coef(mle.weib)["a"]
weib.b <- coef(mle.weib)["b"]

mle.lnorm <-  mle2(NLL.lnorm, start=list(a=2, b=0.5), method="BFGS")
lnorm.a <- coef(mle.lnorm)["a"]
lnorm.b <- coef(mle.lnorm)["b"]

mle.wald <- mle2(NLL.wald, start=list(a=10, b=10), method="BFGS")
wald.a <- coef(mle.wald)["a"]
wald.b <- coef(mle.wald)["b"]




mle.orange1 <- mle.lnorm
################################################################################




pinkMod <- NULL
orangeMod <- NULL
pink <- pink2[,3:4] 
orange <- orange2[,3:4] 

##########################
for (i in 1:nrow(pink)){
  dists <- pink[i,2]
  reps <- pink[i,1]
  pinkVec <- rep(dists, reps)
  pinkMod <- c(pinkMod, pinkVec)
}

dists <- pinkMod

mle.NE <- mle2(NLL.NE, start=list(a=10), method="BFGS")
ne.disp <- coef(mle.NE)["a"]
#
#stud.t <- mle2(NLL.st, start=list(a=10,b=1), method="BFGS")
#st.scale <- coef(stud.t)["a"]
#st.shape <- coef(stud.t)["b"]
  
mle.weib <-  mle2(NLL.weib, start=list(a=10, b=1), method="BFGS")
weib.a <- coef(mle.weib)["a"]
weib.b <- coef(mle.weib)["b"]

mle.lnorm <-  mle2(NLL.lnorm, start=list(a=1, b=1), method="BFGS")
lnorm.a <- coef(mle.lnorm)["a"]
lnorm.b <- coef(mle.lnorm)["b"]

mle.wald <- mle2(NLL.wald, start=list(a=20, b=100), method="BFGS")
wald.a <- coef(mle.wald)["a"]
wald.b <- coef(mle.wald)["b"]

print(AICtab(mle.NE, mle.weib, mle.wald, mle.lnorm,
delta=TRUE, weights=TRUE, sort=TRUE))

mle.pink2 <- mle.wald

##########################
for (i in 1:nrow(orange)){
  dists <- orange[i,2]
  reps <- orange[i,1]
  orangeVec <- rep(dists, reps)
  orangeMod <- c(orangeMod, orangeVec)
}

dists <- orangeMod

mle.NE <- mle2(NLL.NE, start=list(a=10), method="BFGS")
ne.disp <- coef(mle.NE)["a"]

#stud.t <- mle2(NLL.st, start=list(a=10,b=1), method="BFGS")
#st.scale <- coef(stud.t)["a"]
#st.shape <- coef(stud.t)["b"]
  
mle.weib <-  mle2(NLL.weib, start=list(a=10, b=1), method="BFGS")
weib.a <- coef(mle.weib)["a"]
weib.b <- coef(mle.weib)["b"]

mle.lnorm <-  mle2(NLL.lnorm, start=list(a=1, b=1), method="BFGS")
lnorm.a <- coef(mle.lnorm)["a"]
lnorm.b <- coef(mle.lnorm)["b"]

mle.wald <- mle2(NLL.wald, start=list(a=10, b=10), method="BFGS")
wald.a <- coef(mle.wald)["a"]
wald.b <- coef(mle.wald)["b"]

print(AICtab(mle.NE, mle.weib, mle.wald, mle.lnorm,
delta=TRUE, weights=TRUE, sort=TRUE))

mle.orange2 <- mle.wald

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