Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
0 parents
commit e6474c5
Showing
11 changed files
with
536 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,22 @@ | ||
Package: FracKrigingR | ||
Type: Package | ||
Title: Spatial Multivariate Data Modeling | ||
Version: 1.0.0 | ||
Authors@R: c(person("Neringa", "Urbonaite", email = "neringa.urbonaite@mif.vu.lt", role = c("aut", "cre")), | ||
person("Leonidas", "Sakalauskas", email = "leonidas.sakalauskas@mif.vu.lt", role = c("aut"))) | ||
Copyright: Vilnius University Institute of Data Science and Digital | ||
Technologies | ||
Author: Neringa Urbonaite [aut, cre], | ||
Leonidas Sakalauskas [aut] | ||
Maintainer: Neringa Urbonaite <neringa.urbonaite@mif.vu.lt> | ||
Description: Aim is to provide fractional Brownian vector field generation algorithm, Hurst parameter estimation method and fractional kriging model for multivariate data modeling. | ||
License: GPL-2 | ||
Encoding: UTF-8 | ||
URL: https://github.com/NidaGreen/FracKriging | ||
Imports: psych, clusterGeneration, graphics, stats | ||
Suggests: knitr, gstat, sp, rmarkdown, raster | ||
RoxygenNote: 7.1.2 | ||
NeedsCompilation: no | ||
Packaged: 2021-11-05 13:42:32 UTC; nerin | ||
Repository: CRAN | ||
Date/Publication: 2021-11-08 08:40:08 UTC |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,10 @@ | ||
32294b9e9a7f24a92988e7df36ad82f5 *DESCRIPTION | ||
fe09b42dbe196b45447dbdda66b56dba *NAMESPACE | ||
b08b0433c4b7157251526e7c53186fc5 *R/FracField.R | ||
541d30f9ba5af5b93719d46b4ff8569e *R/FracKrig.R | ||
ced8e8af452d155f5223e51372e2a01d *R/FracMatrix.R | ||
fe7b2558ab5f66425c67e2c0037228ed *R/MaxLikelihood.R | ||
44fce02cf14a6696bd669d1356a8e419 *man/FracField.Rd | ||
ccf970a32272d9283dab2747a8056d06 *man/FracKrig.Rd | ||
ab6009b9ab364363bbee028ee6ca657a *man/FracMatrix.Rd | ||
5e5b9f3db7643cedf7d0eaefe6a36ed2 *man/MaxLikelihood.Rd |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,8 @@ | ||
|
||
|
||
export(FracField) | ||
export(FracKrig) | ||
export(FracMatrix) | ||
export(MaxLikelihood) | ||
importFrom("graphics", "mtext", "points") | ||
importFrom("stats", "rnorm") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,60 @@ | ||
#' @title FracField | ||
#' | ||
#' @description Generates fractional Brownian vector field data | ||
#' | ||
#' @param K number of observations | ||
#' @param m number of criteria | ||
#' @param H Hurst parameter (a real in interval [0,1)) | ||
#' @param X Coordinates | ||
|
||
#' @return Returns a fractional Brownian vector field matrix. | ||
#' | ||
#' @examples | ||
#' # Load FracKrigingR library | ||
#' library(FracKrigingR) | ||
#' # generate Coordinates | ||
#' p=2; K=10; | ||
#' X<-matrix(0,ncol=p, nrow=K) | ||
#' for(j in 1:p){ | ||
#' for(i in 1:K){ | ||
#' X[i,j] = rnorm(1, 0, 1) | ||
#' } | ||
#' } | ||
#' # generate fractional Brownian vector field | ||
#' H = 0.5; m = 3 | ||
#' FracField(K,m,H,X) | ||
#' | ||
#' @importFrom graphics mtext | ||
#' | ||
#' @importFrom points stats rnorm | ||
#' @export | ||
|
||
|
||
|
||
|
||
FracField<- function(K,m,H,X){ | ||
suppressWarnings({ | ||
positive_definite_matrix <- (clusterGeneration::genPositiveDefMat(dim = m)$Sigma) | ||
cholecky_matrix<- t(chol(positive_definite_matrix)) | ||
|
||
E<-as.matrix( c(rep(c(1), times = K))) | ||
|
||
gB<-function(H){ | ||
gg<-((E%*%t(E))-((FracMatrix(H,K,X)*c(t(E)%*%solve(FracMatrix(H,K,X))%*%E))/2)) | ||
gg | ||
|
||
} | ||
|
||
vector<- function(m,K){ | ||
ksi<-matrix(0,ncol =m,nrow=K) | ||
for (l in seq(from =1,to =m, by=1)){ksi[,l]<-rnorm(K,0,1)}; | ||
ksi | ||
} | ||
|
||
y<-t(cholecky_matrix%*%t(vector(m,K))%*%((chol(gB(H))))) | ||
|
||
return(y) | ||
}) | ||
} | ||
|
||
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,106 @@ | ||
#' @title FracKrig | ||
#' | ||
#' @description Performs extrapolation for spatial multivariate data | ||
#' | ||
#' @param X Coordinates | ||
#' @param Z observations | ||
#' @param Xnew Coordinates of points where the prognosis should be made | ||
#' @param H Hurst parameter (a real in interval [0,1)) | ||
#' | ||
#' @return Returns a matrix of fractional kriging prognosis. | ||
#' | ||
#' | ||
#' @examples | ||
#' | ||
#' library(sp) | ||
#' library(gstat) | ||
#' data(meuse) | ||
#' xy<-cbind(meuse$x,meuse$y) | ||
|
||
#' X<-xy[1:50,] | ||
#' min_max_norm <- function(x) { | ||
#' (x - min(x)) / (max(x) - min(x)) | ||
#' } | ||
#' normalize <- function(x) { | ||
#' return ((x - min(x)) / (max(x) - min(x))) | ||
#' } | ||
|
||
|
||
#' dat<-cbind(meuse[3],meuse[4],meuse[5]) | ||
#' data<-dat[51:100,] | ||
#' zz1 <- as.data.frame(lapply(dat, normalize)) | ||
#' data1=as.data.frame(lapply(as.data.frame(data), normalize)) | ||
#' Z<-as.matrix(zz1[1:50,]) | ||
|
||
# Load FracKrigingR library | ||
#' library(FracKrigingR) | ||
#' K<-50 | ||
|
||
#' #Hurst parameter estimation | ||
#' H<-0.2 | ||
#' Xnew<-xy[51:100,] | ||
#' results<- FracKrig(X,Z,Xnew,H) | ||
#' denormalize <- function(x, bottom, top){ | ||
#' (top - bottom) * x + bottom | ||
#' } | ||
|
||
#'z1 = denormalize( | ||
#' results[,1], top = max(data[,1]), bottom = min(data[,1]) | ||
#') | ||
#'z2 = denormalize( | ||
#' results[,2], top = max(data[,2]), bottom = min(data[,2]) | ||
#') | ||
#'z3 = denormalize( | ||
#' results[,3], top = max(data[,3]), bottom = min(data[,3]) | ||
#') | ||
|
||
|
||
#'RMSE<-function(z,prognosis){ | ||
#' rmse<-sqrt(((1/(length(z))))*sum((z-prognosis)^2)) | ||
#' rmse | ||
#'} | ||
|
||
#'Cd<-RMSE(data[,1],z1) | ||
#'Cu<-RMSE(data[,2],z2) | ||
#'Pb<-RMSE(data[,3],z3) | ||
|
||
|
||
#'Cd | ||
#'Cu | ||
#'Pb | ||
#' | ||
#' @export | ||
#' | ||
|
||
|
||
FracKrig<- function(X,Z,Xnew,H){ | ||
X<-as.matrix(X) | ||
Z<-as.matrix(Z) | ||
Xnew<-as.matrix(Xnew) | ||
E<-c(rep(1,K)) | ||
E<-as.matrix(E) | ||
K<-nrow(X) | ||
p<-ncol(X) | ||
P<-nrow(Xnew) | ||
m<-ncol(Z) | ||
|
||
prognosis<-function(xpr,H){ | ||
|
||
b<-c() | ||
for (i in seq(from =1, to=K, by=1)){ | ||
b<-c(b,(((xpr - t(X[i,]))%*%t(xpr-t(X[i,])))^H)) | ||
} | ||
|
||
prog<-(as.matrix(t(Z)%*%solve(FracMatrix(H,K,X))))%*%((as.matrix(b))+(E%*%(((1-c(t(E)%*%solve(FracMatrix(H,K,X))%*%(as.matrix(b))))/c(t(E)%*%solve(FracMatrix(H,K,X))%*%E))))) | ||
prog | ||
|
||
} | ||
|
||
|
||
final_result = matrix(c(0), nrow= P, ncol =m) | ||
for (i in seq(from =1, to=P, by=1)){ | ||
|
||
final_result[i,]<-t(prognosis(Xnew[i,],H)) | ||
} | ||
final_result | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,37 @@ | ||
#' @title FracMatrix | ||
#' | ||
#' @description Fractional distance matrix | ||
#' | ||
#' @param H Hurst parameter (a real in interval [0,1)) | ||
#' @param K number of observations | ||
#' @param X Coordinates | ||
#' | ||
#' @return Returns a fractional distance matrix, which depends on the Hurst parameter. | ||
#' | ||
#' @examples | ||
#' # Load FracKrigingR library | ||
#' library(FracKrigingR) | ||
#' #Fractional Brownian vector field | ||
#' K = 10; H = 0.5; p = 2 | ||
#' #Generate coordinates | ||
#' X<-matrix(0,ncol=p, nrow=K) | ||
#' for(j in 1:p){ | ||
#' for(i in 1:K){ | ||
#' X[i,j] = rnorm(1, 0, 1) | ||
#' } | ||
#' } | ||
#' FracMatrix(H, K, X) | ||
#' @export | ||
|
||
|
||
|
||
|
||
FracMatrix<-function(H, K, X){ | ||
a<-matrix(0,nrow=K,ncol=K) | ||
for (i in 1:K){ | ||
for (j in 1:K){ | ||
a[i,j]<- ((t(X[i,]-X[j,]))%*%((X[i,])-X[j,]))^H | ||
} | ||
} | ||
a | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,110 @@ | ||
#' @title MaxLikelihood | ||
#' | ||
#' @description Maximum likelihood method for Hurst parameter estimation of multivariate data | ||
#' | ||
#' @param X Coordinates | ||
#' @param Z Observations | ||
#' | ||
#' @return Returns the estimate of the Hurst parameter (a real in [0,1)) | ||
#' and a graph indicating the minimized maximum likelihood function with the Hurst parameter. | ||
#' @examples | ||
#' # Load FracKrigingR library | ||
#' library(FracKrigingR) | ||
#' # generate Coordinates | ||
#' p<-2; K<-20; | ||
#' X<-matrix(0,ncol=p, nrow=K) | ||
#' for(j in 1:p){ | ||
#' for(i in 1:K){ | ||
#' X[i,j] = rnorm(1, 0, 1) | ||
#' } | ||
#' } | ||
#' # generate fractional Brownian vector field | ||
#' H <- 0.8; m <- 3 | ||
#' Z<-FracField(K,m,H,X) | ||
#' # Hurst parameter estimation | ||
#' MaxLikelihood(X,Z) | ||
#' | ||
#' @export | ||
|
||
MaxLikelihood<- function(X,Z){ | ||
Z<-as.matrix(Z) | ||
K = nrow(X) | ||
m = ncol(Z) | ||
E<-c(rep(1,K)) | ||
E<-as.matrix(E) | ||
|
||
|
||
golden.section.search = function(f, lower.bound, upper.bound, tolerance) | ||
{ | ||
golden.ratio = 2/(sqrt(5) + 1) | ||
|
||
x1 = upper.bound - golden.ratio*(upper.bound - lower.bound) | ||
x2 = lower.bound + golden.ratio*(upper.bound - lower.bound) | ||
|
||
f1 = f(x1) | ||
f2 = f(x2) | ||
|
||
ite = 0 | ||
while (abs(upper.bound - lower.bound) > tolerance) | ||
{ | ||
ite = ite + 1 | ||
if (f2 > f1) | ||
{ | ||
upper.bound = x2 | ||
x2 = x1 | ||
f2 = f1 | ||
x1 = upper.bound - golden.ratio*(upper.bound - lower.bound) | ||
f1 = f(x1) | ||
} | ||
else | ||
{ | ||
lower.bound = x1 | ||
x1 = x2 | ||
f1 = f2 | ||
x2 = lower.bound + golden.ratio*(upper.bound - lower.bound) | ||
f2 = f(x2) | ||
} | ||
} | ||
minn = (lower.bound + upper.bound)/2 | ||
minn | ||
} | ||
|
||
|
||
|
||
|
||
A<-function(H){ | ||
|
||
a<-matrix(0,nrow=K,ncol=K) | ||
for (i in 1:K){ | ||
for (j in 1:K){ | ||
a[i,j]<- ((t(X[i,]-X[j,]))%*%((X[i,])-X[j,]))^H | ||
} | ||
} | ||
a | ||
} | ||
|
||
bc<-function(H){ | ||
bc<-((1/K)*(((((t(Z)%*%solve(A(H))%*%E)%*%t((t(Z)%*%solve(A(H))%*%E)))/c(t(E)%*%solve(A(H))%*%E)))-(t(Z)%*%solve(A(H))%*%Z))) | ||
bc | ||
} | ||
|
||
Ff<-function(H){ | ||
f<-(log(det(bc(H)))/m)+log(abs(det(A(H))))/K | ||
f | ||
} | ||
|
||
|
||
t<-seq(from=0.01,to=0.9,by=0.01) | ||
vect<-c() | ||
for (i in t){ | ||
vect<-c(vect,Ff(i)) | ||
} | ||
gsc<-golden.section.search(Ff, 0.01, 0.9, 1e-5) | ||
ff<-Ff(gsc) | ||
|
||
plot(t,vect,type="l", col="darkblue", lwd=1, xlab="", ylab="",cex.axis=1.3) | ||
points(gsc,ff,pch=8,col="darkgreen",cex=2) | ||
mtext(text = "Hurst Parameter", side = 1,line = 2.4,font=1,cex=1.4) | ||
mtext(text = "Maximum Likelihood", side = 2, line = 2.4,font=1,cex=1.4) | ||
gsc | ||
} |
Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.
Oops, something went wrong.
Oops, something went wrong.