bayes6<-function(n,dataset, U,  sigmastart,xistart,alphastart,  erreta,errxi,erralpha,  sdeta,sdxi)
{

#sigmastart takes two values, as does xistart; erreta takes two values, as does errxi; as do both of sdeta and sdxi

#Length of the dataset (NB: dataset is a 2-column matrix)
k<-length(dataset[,1])

#Vectors to take candidate values:
caneta_x = vector("numeric")  #eta = log(sigma)
caneta_y = vector("numeric")
canxi_x = vector("numeric")
canxi_y = vector("numeric")
canalpha = vector("numeric")

#Initialise the candidate values at their starting values:
caneta_x[1] = log(sigmastart[1])
caneta_y[1] = log(sigmastart[2])
canxi_x[1] = xistart[1]
canxi_y[1] = xistart[2]
canalpha[1] = alphastart

#Vectors to take final values:
x_eta = vector("numeric")
y_eta = vector("numeric")
x_xi = vector("numeric")
y_xi = vector("numeric")
z = vector("numeric")

#Vectors to take the acceptance probabilities:
aprobeta_x = vector("numeric")
aprobeta_y = vector("numeric")
aprobxi_x = vector("numeric")
aprobxi_y = vector("numeric")
aprobalpha = vector("numeric")

#Initialise the final values vectors at the first candidates:
x_eta[1] = caneta_x[1]
y_eta[1] = caneta_y[1]
x_xi[1] = canxi_x[1]
y_xi[1] = canxi_y[1]
z[1] = canalpha[1]

#Now the function which returns the log-likelihood for (log(sigma), xi, alpha):
loglik = function(data, thresh, ETA_x, ETA_y, XI_x, XI_y, ALPHA)
{
        u_x = thresh[1]
        u_y = thresh[2]
 
        if (ALPHA<0.00001 | ALPHA>0.99999)return(-1000000)
        n = dim(data)[1]        
        lambda_x = sum(data[,1] > u_x)/n
        lambda_y = sum(data[,2] > u_y)/n
        quad1 = data[data[,1] > u_x & data[,2] > u_y,]
        quad2 = data[data[,1] > u_x & data[,2] <= u_y,]
        quad3 = data[data[,1] <= u_x & data[,2] > u_y,]
	#Add a bit in for "quad 4" (both below their thresholds)
	quad4 = data[data[,1] <= u_x & data[,2] <= u_y,] 
	#Add in a trap just in case there's only one observation in any of the quadrants - in which case the routine fails as it doesn't see "quad" as a matrix, and so can't read the columns...

	#if(!is.matrix(quad1)){quad1 = t(as.matrix(quad1))}
	#if(!is.matrix(quad2)){quad2 = t(as.matrix(quad2))}
	#if(!is.matrix(quad3)){quad3 = t(as.matrix(quad3))}
	#if(!is.matrix(quad4)){quad4 = t(as.matrix(quad4))}

        core1x = 1 + XI_x*(quad1[,1]-u_x)/exp(ETA_x)
        core1y = 1 + XI_y*(quad1[,2]-u_y)/exp(ETA_y)
        core2x = 1 + XI_x*(quad2[,1]-u_x)/exp(ETA_x)
        core3y = 1 + XI_y*(quad3[,2]-u_y)/exp(ETA_y)
        if(min(core1x)<0.00001)return(-1000000)
        if(min(core1y)<0.00001)return(-1000000)
        if(min(core2x)<0.00001)return(-1000000)
        if(min(core3y)<0.00001)return(-1000000)

        xtil1 = -(log(1-lambda_x * (1+XI_x*(quad1[,1]-u_x)/exp(ETA_x))**(-1/XI_x)))**(-1)
        ytil1 = -(log(1-lambda_y * (1+XI_y*(quad1[,2]-u_y)/exp(ETA_y))**(-1/XI_y)))**(-1)
        dxtil1bydx = ((lambda_x/exp(ETA_x))*core1x**((-1/XI_x)-1))*((log(1-lambda_x*core1x**(-1/XI_x)))**(-2))/(1-lambda_x*core1x**(-1/XI_x)) 
        dytil1bydy = ((lambda_y/exp(ETA_y))*core1y**((-1/XI_y)-1))*((log(1-lambda_y*core1y**(-1/XI_y)))**(-2))/(1-lambda_y*core1y**(-1/XI_y))
        bit1 = xtil1**(-1/ALPHA) + ytil1**(-1/ALPHA)
        loglik1 = -bit1**ALPHA + ((-1/ALPHA)-1)*(log(xtil1) + log(ytil1)) + (ALPHA-2)*log(bit1) +log(bit1**ALPHA -1 +1/ALPHA) + log(dxtil1bydx) + log(dytil1bydy)
                
        xtil2 = -log((1-lambda_x * (1+XI_x*(quad2[,1]-u_x)/exp(ETA_x))**(-1/XI_x)))**(-1)
        ytil2 = -log(1-lambda_y)**(-1)
        dxtil2bydx = ((lambda_x/exp(ETA_x))*core2x**((-1/XI_x)-1))*((log(1-lambda_x*core2x**(-1/XI_x)))**(-2))/(1-lambda_x*core2x**(-1/XI_x)) 
        bit2 = xtil2**(-1/ALPHA) + ytil2**(-1/ALPHA)
        loglik2 = -bit2**ALPHA + ((-1/ALPHA) - 1)*log(xtil2) + (ALPHA-1)*log(bit2) + log(dxtil2bydx)
        
        xtil3 = -log(1-lambda_x)**(-1)
        ytil3 = -log((1-lambda_y * (1+XI_y*(quad3[,2]-u_y)/exp(ETA_y))**(-1/XI_y)))**(-1)
        dytil3bydy = ((lambda_y/exp(ETA_y))*core3y**((-1/XI_y)-1))*((log(1-lambda_y*core3y**(-1/XI_y)))**(-2))/(1-lambda_y*core3y**(-1/XI_y)) 
        bit3 = xtil3**(-1/ALPHA) + ytil3**(-1/ALPHA)
        loglik3 = - bit3**ALPHA + ((-1/ALPHA) - 1)*log(ytil3) + (ALPHA-1)*log(bit3) + log(dytil3bydy)

	xtil4 = -log(1-lambda_x)**(-1)
	ytil4 = -log(1-lambda_y)**(-1)
	loglik4 = length(quad4[,1])*exp(-(((xtil4^(-1/ALPHA))+(ytil4^(-1/ALPHA)))^(ALPHA)))

        llm = (sum(loglik1)+sum(loglik2)+sum(loglik3)+loglik4)
        return(llm)
}


for(i in 2:n){
#print(i)
caneta_x[i]<-x_eta[i-1]+rnorm(1,0,erreta[1])

#(data,thresh, ETA_x, ETA_y, XI_x, XI_y, ALPHA)
#First update done.  Check!!!!!

likely<-exp((loglik(dataset, U, caneta_x[i],y_eta[i-1],x_xi[i-1],y_xi[i-1],z[i-1]))-(loglik(dataset,U,x_eta[i-1],y_eta[i-1],x_xi[i-1],y_xi[i-1],z[i-1])))

aprobeta_x[i]<-min(1,(likely*((dnorm(caneta_x[i],0,sdeta[1]))*(dnorm(y_eta[i-1],0,sdeta[2]))* (dnorm(x_xi[i-1],0,sdxi[1])) *(dnorm(y_xi[i-1],0,sdxi[2]))*(dunif(z[i-1])) )/((dnorm(x_eta[i-1],0,sdeta[1]))*(dnorm(y_eta[i-1],0,sdeta[2]))* (dnorm(x_xi[i-1],0,sdxi[1])) *(dnorm(y_xi[i-1],0,sdxi[2]))*(dunif(z[i-1])) )))

print(c(i, caneta_x[i], x_eta[i-1], y_eta[i-1], x_xi[i-1],y_xi[i-1], z[i-1]))

u<-runif(1)
if(u<aprobeta_x[i]){x_eta[i]<-caneta_x[i]}
if(u>=aprobeta_x[i]){x_eta[i]<-x_eta[i-1]}
#####################################################################
caneta_y[i]<-y_eta[i-1]+rnorm(1,0,erreta[2])
likely2<-exp((loglik(dataset, U, x_eta[i],caneta_y[i],x_xi[i-1],y_xi[i-1],z[i-1]))-(loglik(dataset,U,x_eta[i],y_eta[i-1],x_xi[i-1],y_xi[i-1],z[i-1])))

aprobeta_y[i]<-min(1,(likely2*((dnorm(x_eta[i],0,sdeta[1]))*(dnorm(caneta_y[i],0,sdeta[2]))* (dnorm(x_xi[i-1],0,sdxi[1])) *(dnorm(y_xi[i-1],0,sdxi[2]))*(dunif(z[i-1])) )/((dnorm(x_eta[i],0,sdeta[1]))*(dnorm(y_eta[i-1],0,sdeta[2]))* (dnorm(x_xi[i-1],0,sdxi[1])) *(dnorm(y_xi[i-1],0,sdxi[2]))*(dunif(z[i-1])) )))

print(c(i, caneta_y[i], x_eta[i], y_eta[i-1], x_xi[i-1],y_xi[i-1], z[i-1]))

u<-runif(1)
if(u<aprobeta_y[i]){y_eta[i]<-caneta_y[i]}
if(u>=aprobeta_y[i]){y_eta[i]<-y_eta[i-1]}

#########################################################################
canxi_x[i]<-x_xi[i-1]+rnorm(1,0,errxi[1])

likely3<-exp((loglik(dataset, U, x_eta[i],y_eta[i],canxi_x[i],y_xi[i-1],z[i-1]))-(loglik(dataset,U,x_eta[i],y_eta[i],x_xi[i-1],y_xi[i-1],z[i-1])))

aprobxi_x[i]<-min(1,(likely3*((dnorm(x_eta[i],0,sdeta[1]))*(dnorm(y_eta[i],0,sdeta[2]))* (dnorm(canxi_x[i],0,sdxi[1])) *(dnorm(y_xi[i-1],0,sdxi[2]))*(dunif(z[i-1])) )/((dnorm(x_eta[i],0,sdeta[1]))*(dnorm(y_eta[i],0,sdeta[2]))* (dnorm(x_xi[i-1],0,sdxi[1])) *(dnorm(y_xi[i-1],0,sdxi[2]))*(dunif(z[i-1])) )))

print(c(i, canxi_x[i], x_eta[i], y_eta[i], x_xi[i-1],y_xi[i-1], z[i-1]))

u<-runif(1)
if(u<aprobxi_x[i]){x_xi[i]<-canxi_x[i]}
if(u>=aprobxi_x[i]){x_xi[i]<-x_xi[i-1]}

#############################################################################
canxi_y[i]<-y_xi[i-1]+rnorm(1,0,errxi[2])

likely4<-exp((loglik(dataset, U, x_eta[i],y_eta[i],x_xi[i],canxi_y[i],z[i-1]))-(loglik(dataset,U,x_eta[i],y_eta[i],x_xi[i],y_xi[i-1],z[i-1])))

aprobxi_y[i]<-min(1,(likely4*((dnorm(x_eta[i],0,sdeta[1]))*(dnorm(y_eta[i],0,sdeta[2]))* (dnorm(x_xi[i],0,sdxi[1])) *(dnorm(canxi_y[i],0,sdxi[2]))*(dunif(z[i-1])) )/((dnorm(x_eta[i],0,sdeta[1]))*(dnorm(y_eta[i],0,sdeta[2]))* (dnorm(x_xi[i],0,sdxi[1])) *(dnorm(y_xi[i-1],0,sdxi[2]))*(dunif(z[i-1])) )))

print(c(i, canxi_y[i], x_eta[i], y_eta[i], x_xi[i],y_xi[i-1], z[i-1]))

u<-runif(1)
if(u<aprobxi_y[i]){y_xi[i]<-canxi_y[i]}
if(u>=aprobxi_y[i]){y_xi[i]<-y_xi[i-1]}

##########################################################################
canalpha[i]<-z[i-1]+rnorm(1,0,erralpha)

likely5<-exp((loglik(dataset, U, x_eta[i],y_eta[i],x_xi[i],y_xi[i],canalpha[i]))-(loglik(dataset,U,x_eta[i],y_eta[i],x_xi[i],y_xi[i],z[i-1])))

aprobalpha[i]<-min(1,(likely5*((dnorm(x_eta[i],0,sdeta[1]))*(dnorm(y_eta[i],0,sdeta[2]))* (dnorm(x_xi[i],0,sdxi[1])) *(dnorm(y_xi[i],0,sdxi[2]))*(dunif(canalpha[i])) )/((dnorm(x_eta[i],0,sdeta[1]))*(dnorm(y_eta[i],0,sdeta[2]))* (dnorm(x_xi[i],0,sdxi[1])) *(dnorm(y_xi[i],0,sdxi[2]))*(dunif(z[i-1])) )))

print(c(i, canalpha[i], x_eta[i], y_eta[i], x_xi[i],y_xi[i], z[i-1]))

u<-runif(1)
if(u<aprobalpha[i]){z[i]<-canalpha[i]}
if(u>=aprobalpha[i]){z[i]<-z[i-1]}
#print(i)

}


#aprobeta_x<-aprobeta_x[!is.na(aprobeta_x)]
#aprobeta_y<-aprobeta_y[!is.na(aprobeta_y)]
#aprobxi_x<-aprobxi_x[!is.na(aprobxi_x)]
#aprobxi_y<-aprobxi_y[!is.na(aprobxi_y)]
#aprobalpha<-aprobalpha[!is.na(aprobalpha)]

results=matrix(ncol=10,nrow=n)
results[,1]=x_eta
results[,2]=y_eta
results[,3]=x_xi
results[,4]=y_xi
results[,5]=z
results[,6]=aprobeta_x
results[,7]=aprobeta_y
results[,8]=aprobxi_x
results[,9]=aprobxi_y
results[,10]=aprobalpha
return(results)
}

##TEST:
library(ismev)
library(evd)

data(wavesurge)
attach(wavesurge)
test = wavesurge

thresh.p = 0.85

#test = rbvevd(1000, dep=0.5, model="log")

U = as.numeric(c(quantile(test[,1],thresh.p), quantile(test[,2],thresh.p)))

SIGMASTART = c(sd(test[,1]), sd(test[,2]))
XISTART = c(0.1,0.1)
ALPHASTART = 0.5
ERRETA = c(0.2,0.15)
ERRXI = c(0.1, 0.1)
ERRALPHA = 0.1
SDETA = c(1000,1000)
SDXI = c(100, 100)

lee = bayes6(1000, test, U, SIGMASTART, XISTART, ALPHASTART,  ERRETA, ERRXI, ERRALPHA, SDETA, SDXI)


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