[aggexp.git] / pkg / R / z_plotHelper.R

#' @include z_runAlgorithm.R

#' @title Get best expert index
#'
#' @description Return the weights corresponding to the best expert (...0,1,0...)
#'
#' @param r Output of \code{\link{runAlgorithm}}
#'
#' @export
getBestExpert = function(r)
{
	X = as.matrix(r$data[,names(r$data) %in% r$experts])
	Y = r$data[,"Measure"]

	bestIndex = which.min(colMeans(abs(X - Y)^2, na.rm=TRUE))
	res = rep(0.0, length(r$experts))
	res[bestIndex] = 1.0
	return (res)
}

#' @title Get best convex combination
#'
#' @description Return the weights p minimizing the quadratic error ||X*p-Y||^2 under convexity contraint.
#'
#' @param r Output of \code{\link{runAlgorithm}}
#'
#' @export
getBestConvexCombination = function(r)
{
	X = r$data[,r$experts]
	Y = as.double(r$data[,"Measure"])
	indices = getNoNAindices(X) & !is.na(Y)
	X = as.matrix(X[indices,])
	Y = Y[indices]

	K = length(r$experts)
	return (constrOptim(theta=rep(1.0/K,K),
		method="Nelder-Mead", #TODO: others not better... why?
		f=function(p){return(sum((X%*%p-Y)^2))}, 
		grad=NULL, #function(p){return(2.*t(X)%*%(X%*%p-Y))}, 
		ui=rbind(rep(1.,K),rep(-1.,K),diag(K)), ci=c(0.99999,-1.00001, rep(0.,K)), 
		control=list(ndeps=1e-3,maxit=10000))$par)
}

#' @title Get best linear combination
#'
#' @description Return the weights u minimizing the quadratic error ||r$X*u-r$Y||^2
#'
#' @param r Output of \code{\link{runAlgorithm}}
#'
#' @export
getBestLinearCombination = function(r)
{
	X = r$data[,r$experts]
	Y = r$data[,"Measure"]
	indices = getNoNAindices(X) & !is.na(Y)
	X = as.matrix(X[indices,])
	Y = Y[indices]

	return (mpPsInv(X) %*% Y)
}

#' @title Get statistical indicators
#'
#' @description Return respectively the TS, FA, MA, RMSE, EV indicators in a list.
#'
#' @param r Output of \code{\link{runAlgorithm}}
#' @param thresh Threshold to compute alerts indicators.
#' @param station Name or index of the station to consider. Default: the first one
#' @param noNA TRUE to show only errors associated with full lines (old behavior)
#'
#' @export
getIndicators = function(r, thresh, station=1, noNA=TRUE)
{
	if (is.character(station))
		station = match(station, r$stations)

	#TODO: duplicated block (same in plotCloud())
	XY = subset(r$data, subset = (Station == station), select = c(r$experts,"Measure","Prediction"))
	Y = XY[,"Measure"]
	hatY = XY[,"Prediction"]
	indices = !is.na(Y) & !is.na(hatY)
	if (noNA)
	{
		X = XY[,names(XY) %in% r$experts]
		indices = indices & getNoNAindices(X)
	}
	Y = Y[indices]
	hatY = hatY[indices]

	RMSE = round(sqrt(sum((Y - hatY)^2) / length(Y)),2)
	EV = round(1 - var(Y-hatY) / var(Y), 2)
	A = sum(hatY >= thresh & Y >= thresh, na.rm=TRUE) #right alarm
	B = sum(hatY >= thresh & Y < thresh, na.rm=TRUE) #false alarm
	C = sum(hatY < thresh & Y >= thresh, na.rm=TRUE) #missed alert
	TS = round(A/(A+B+C),2)
	FA = B/(A+B)
	MA = C/(A+C)
	return (list("TS"=TS, "FA"=FA, "MA"=MA, "RMSE"=RMSE, "EV"=EV))
}
Commit	Line	Data
a961f8a1 BA	1	#' @include z_runAlgorithm.R
	2
	3	#' @title Get best expert index
	4	#'
	5	#' @description Return the weights corresponding to the best expert (...0,1,0...)
	6	#'
	7	#' @param r Output of \code{\link{runAlgorithm}}
	8	#'
	9	#' @export
	10	getBestExpert = function(r)
	11	{
	12	X = as.matrix(r$data[,names(r$data) %in% r$experts])
	13	Y = r$data[,"Measure"]
	14
	15	bestIndex = which.min(colMeans(abs(X - Y)^2, na.rm=TRUE))
	16	res = rep(0.0, length(r$experts))
	17	res[bestIndex] = 1.0
	18	return (res)
	19	}
	20
	21	#' @title Get best convex combination
	22	#'
	23	#' @description Return the weights p minimizing the quadratic error \|\|X*p-Y\|\|^2 under convexity contraint.
	24	#'
	25	#' @param r Output of \code{\link{runAlgorithm}}
	26	#'
	27	#' @export
	28	getBestConvexCombination = function(r)
	29	{
	30	X = r$data[,r$experts]
	31	Y = as.double(r$data[,"Measure"])
	32	indices = getNoNAindices(X) & !is.na(Y)
	33	X = as.matrix(X[indices,])
	34	Y = Y[indices]
	35
	36	K = length(r$experts)
	37	return (constrOptim(theta=rep(1.0/K,K),
	38	method="Nelder-Mead", #TODO: others not better... why?
	39	f=function(p){return(sum((X%*%p-Y)^2))},
	40	grad=NULL, #function(p){return(2.t(X)%%(X%*%p-Y))},
	41	ui=rbind(rep(1.,K),rep(-1.,K),diag(K)), ci=c(0.99999,-1.00001, rep(0.,K)),
	42	control=list(ndeps=1e-3,maxit=10000))$par)
	43	}
	44
	45	#' @title Get best linear combination
	46	#'
	47	#' @description Return the weights u minimizing the quadratic error \|\|r$X*u-r$Y\|\|^2
	48	#'
	49	#' @param r Output of \code{\link{runAlgorithm}}
	50	#'
	51	#' @export
	52	getBestLinearCombination = function(r)
	53	{
	54	X = r$data[,r$experts]
	55	Y = r$data[,"Measure"]
	56	indices = getNoNAindices(X) & !is.na(Y)
	57	X = as.matrix(X[indices,])
	58	Y = Y[indices]
	59
	60	return (mpPsInv(X) %*% Y)
	61	}
	62
	63	#' @title Get statistical indicators
	64	#'
65	#' @description Return respectively the TS, FA, MA, RMSE, EV indicators in a list.
66	#'
67	#' @param r Output of \code{\link{runAlgorithm}}
68	#' @param thresh Threshold to compute alerts indicators.
69	#' @param station Name or index of the station to consider. Default: the first one
70	#' @param noNA TRUE to show only errors associated with full lines (old behavior)
71	#'
72	#' @export
73	getIndicators = function(r, thresh, station=1, noNA=TRUE)
74	{
75	if (is.character(station))
76	station = match(station, r$stations)
77
78	#TODO: duplicated block (same in plotCloud())
79	XY = subset(r$data, subset = (Station == station), select = c(r$experts,"Measure","Prediction"))
80	Y = XY[,"Measure"]
81	hatY = XY[,"Prediction"]
82	indices = !is.na(Y) & !is.na(hatY)
83	if (noNA)
84	{
85	X = XY[,names(XY) %in% r$experts]
86	indices = indices & getNoNAindices(X)
87	}
88	Y = Y[indices]
89	hatY = hatY[indices]
90
91	RMSE = round(sqrt(sum((Y - hatY)^2) / length(Y)),2)
92	EV = round(1 - var(Y-hatY) / var(Y), 2)
93	A = sum(hatY >= thresh & Y >= thresh, na.rm=TRUE) #right alarm
94	B = sum(hatY >= thresh & Y < thresh, na.rm=TRUE) #false alarm
95	C = sum(hatY < thresh & Y >= thresh, na.rm=TRUE) #missed alert
96	TS = round(A/(A+B+C),2)
97	FA = B/(A+B)
98	MA = C/(A+C)
99	return (list("TS"=TS, "FA"=FA, "MA"=MA, "RMSE"=RMSE, "EV"=EV))
100	}