DiceKriging/R/logLikFun.R at 8ccc5ef4f187e9cbcf7082344a2868973d5adb31 · IRSN/DiceKriging · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
`logLikFun` <-
function(param, model, envir=NULL) {
  if (model@known.param=="Trend") {
    beta <- model@trend.coef
  } else {
    beta <- NULL
  }
	if (identical(model@case, "LLconcentration_beta_sigma2")) {

		model@covariance <- vect2covparam(model@covariance, param)
		model@covariance@sd2 <- 1		# to get the correlation matrix

		aux <- covMatrix(model@covariance, model@X)

		R <- aux[[1]]
		T <- chol(R)

   	x <- backsolve(t(T), model@y, upper.tri = FALSE)
   	M <- backsolve(t(T), model@F, upper.tri = FALSE)
		z <- compute.z(x=x, M=M, beta=beta)
		sigma2.hat <- compute.sigma2.hat(z)
		logLik <- -0.5*(model@n * log(2*pi*sigma2.hat) + 2*sum(log(diag(T))) + model@n)

		if (!is.null(envir)) {
			envir$T <- T
      envir$R <- R
      envir$z <- z
      envir$sigma2.hat <- sigma2.hat
		}

	} else if (identical(model@case, "LLconcentration_beta")) {
		nparam <- length(param)
                if (class(model@covariance) != "covAdditive0") {
		  model@covariance <- vect2covparam(model@covariance, param[1:(nparam-1)])
		  model@covariance@sd2 <- param[nparam]
    } else {
       model@covariance <- vect2covparam(model@covariance, param)
    }

# 		if (model@covariance@nugget.estim) {
#       coef(model@covariance, type="all") <- param
# 		} else {
#       coef(model@covariance, type="all-nugget") <- param
# 		}
# 		# pb ici : reconnaitre les variables a estimer ?
    aux <- covMatrix(model@covariance, model@X, noise.var=model@noise.var)

		C <- aux[[1]]
		vn <- aux[[2]]
    T <- chol(C)
    x <- backsolve(t(T), model@y, upper.tri = FALSE)
    M <- backsolve(t(T), model@F, upper.tri = FALSE)
		z <- compute.z(x=x, M=M, beta=beta)

    logLik <-  -0.5*(model@n * log(2*pi) + 2*sum(log(diag(T))) + t(z)%*%z)

		if (!is.null(envir)) {
      envir$T <- T
      envir$C <- C
      envir$vn <- vn
      envir$z <- z
		}

	} else if (identical(model@case, "LLconcentration_beta_v_alpha")) {

		nparam <- length(param)

		model@covariance <- vect2covparam(model@covariance, param[1:(nparam-1)])
		model@covariance@sd2 <- 1
		model@covariance@nugget <- 0
		alpha <- param[nparam]

		aux <- covMatrix(model@covariance, model@X)
		R0 <- aux[[1]]-diag(aux[[2]])
		R <- alpha*R0 + (1-alpha)*diag(model@n)

    T <- chol(R)
   	x <- backsolve(t(T), model@y, upper.tri = FALSE)
   	M <- backsolve(t(T), model@F, upper.tri = FALSE)
		z <- compute.z(x=x, M=M, beta=beta)
		v <- compute.sigma2.hat(z)

		logLik <- -0.5*(model@n * log(2*pi*v) + 2*sum(log(diag(T))) + model@n)

		if (!is.null(envir)) {
		  envir$T <- T
		  envir$R0 <- R0
		  envir$v <- v
		  envir$z <- z
		}

	}


	if (model@method=="PMLE") {
			fun <- match.fun(model@penalty$fun)
			penalty <- - model@n * sum(fun(1/model@covariance@range.val^2, model@penalty$value))
			logLik <- logLik + penalty
	}

	return(logLik)

}