Nothing
#' @title Confidence intervals for the Model Parameters
#'
#' @description Computes confidence intervals for the model parameters.
#'
#' @param object a fitted model object of class \code{slm}.
#' @param parm a vector of integer. Specifies the coordinates of the vector of parameters for which a confidence
#' interval will be given. If missing, all parameters are considered.
#' @param level a number between 0 and 1; the confidence level required.
#' @param ... additional argument(s) for methods.
#'
#' @return This function returns the confidence intervals for the parameters of the model.
#'
#' @importFrom stats lm confint qnorm
#'
#' @references
#' E. Caron, J. Dedecker and B. Michel (2019). Linear regression with stationary errors: the R package slm. \emph{arXiv preprint arXiv:1906.06583}.
#' \url{https://arxiv.org/abs/1906.06583}.
#'
#' @seealso
#' \code{\link[stats:confint.lm]{confint.lm}}.
#'
#' @export
#'
#' @examples
#' data("shan")
#' reg1 = slm(shan$PM_Xuhui ~ . , data = shan, method_cov_st = "fitAR", model_selec = -1)
#' confint(reg1, level = 0.8)
#'
#' data("co2")
#' y = as.vector(co2)
#' x = as.vector(time(co2)) - 1958
#' reg2 = slm(y ~ x + I(x^2) + I(x^3) + sin(2*pi*x) + cos(2*pi*x) + sin(4*pi*x) +
#' cos(4*pi*x) + sin(6*pi*x) + cos(6*pi*x) + sin(8*pi*x) + cos(8*pi*x),
#' method_cov_st = "fitAR", model_selec = -1, plot = TRUE)
#' confint(reg2, level = 0.9)
confint.slm = function(object, parm=NULL, level=0.95, ...) {
mod = object$model
names(mod)[1] = "Y"
reglm = lm(Y~., data=mod)
norm_matrix = diag(object@norm_matrix)
coefs = object$coefficients
Cn = diag(cov_matrix_estimator(object))
levels = 1-level
if (is.null(parm)) {
interval = confint(reglm, level=level, ...)
interval[1:dim(interval)[1],1] = coefs - (qnorm(1-levels/2)*sqrt(Cn))/norm_matrix
interval[1:dim(interval)[1],2] = coefs + (qnorm(1-levels/2)*sqrt(Cn))/norm_matrix
} else {
interval = confint(reglm, parm=parm, level=level, ...)
interval[1:dim(interval)[1],1] = coefs[parm] - (qnorm(1-levels/2)*sqrt(Cn[parm]))/norm_matrix[parm]
interval[1:dim(interval)[1],2] = coefs[parm] + (qnorm(1-levels/2)*sqrt(Cn[parm]))/norm_matrix[parm]
}
return(interval)
}
# confint.slm = function(object, parm=NULL, level=0.95, ...) {
# Y = as.matrix(object$model[1])
# design = cbind(rep(1,length(Y)),as.matrix(object$model[-1])) #intercept added in the design X.
# num_var = dim(design)[2]
# norm_matrix = diag(object@norm_matrix)
# coefs = object$coefficients
# Cn = diag(cov_matrix_estimator(object))
# levels = 1-level
# if (is.null(parm)) {
# interval = matrix(0, nrow=num_var, ncol=2, dimnames = list(c("(Intercept)",names(object$model[,-1])),c(levels/2, 1-levels/2)))
# interval[1:num_var,1] = coefs - (qnorm(1-levels/2)*sqrt(Cn))/norm_matrix
# interval[1:num_var,2] = coefs + (qnorm(1-levels/2)*sqrt(Cn))/norm_matrix
# } else {
# interval = matrix(0, nrow=length(parm), ncol=2, dimnames = list(parm,c(levels/2, 1-levels/2)))
# interval[1:length(parm),1] = coefs[parm] - (qnorm(1-levels/2)*sqrt(Cn[parm]))/norm_matrix[parm]
# interval[1:length(parm),2] = coefs[parm] + (qnorm(1-levels/2)*sqrt(Cn[parm]))/norm_matrix[parm]
# }
# return(interval)
# }
#' @title Summarizing Stationary Linear Model Fits
#'
#' @description Summary method for class "\code{slm}".
#'
#' @param object an object of class "\code{slm}", usually, a result of a call to \code{slm}.
#' @param correlation logical. If \code{TRUE}, the correlation matrix of the estimated parameters is returned and printed.
#' @param symbolic.cor logical. If \code{TRUE}, print the correlations in a symbolic form (see \code{\link[stats:symnum]{symnum}}) rather than as numbers.
#' @param ... further arguments passed to or from other methods.
#'
#' @return The function \code{summary.slm} computes and returns a list of summary statistics of the fitted linear model
#' given in \code{object}, using the components (list elements) "\code{call}" and "\code{terms}" from its argument, plus:
#' \item{residuals }{ the residuals, that is response minus fitted values.}
#' \item{coefficients }{ a \eqn{p*4} matrix with columns for the estimated coefficient, its standard error, z-statistic and corresponding (two-sided) p-value.
#' Aliased coefficients are omitted.}
#' \item{aliased }{ named logical vector showing if the original coefficients are aliased.}
#' \item{sigma }{ the square root of the estimated variance of the error process.}
#' \item{df }{ degrees of freedom, a 3-vector \eqn{(p, n-p, p*)}, the first being the number of non-aliased coefficients, the last being the total number of coefficients.}
#' \item{chi2statistic }{ a 2-vector with the value of the chi2-statistic with its degree of freedom.}
#' \item{r.squared }{ \eqn{R^2}, the 'fraction of variance explained by the model'.}
#' \item{cov.unscaled }{ the matrix \eqn{(X^{t} X)^{-1}}.}
#' \item{correlation }{ the correlation matrix corresponding to the above \code{cov.unscaled}, if \code{correlation = TRUE} is specified.}
#' \item{symbolic.cor }{ (only if \code{correlation} is true.) The value of the argument \code{symbolic.cor}.}
#'
#' @import stats
#'
#' @references
#' E. Caron, J. Dedecker and B. Michel (2019). Linear regression with stationary errors: the R package slm. \emph{arXiv preprint arXiv:1906.06583}.
#' \url{https://arxiv.org/abs/1906.06583}.
#'
#' @seealso
#' The model fitting function \code{\link[slm:slm]{slm}}, \code{\link[base:summary]{summary}}.
#'
#' The function \code{\link[stats:coef]{coef}} extracts the matrix of coefficients with standard errors, z-statistics and p-values.
#'
#' @export
#'
#' @examples
#' data("shan")
#' reg1 = slm(shan$PM_Xuhui ~ . , data = shan, method_cov_st = "fitAR", model_selec = -1)
#' summary(reg1)
#'
#' data("co2")
#' y = as.vector(co2)
#' x = as.vector(time(co2)) - 1958
#' reg2 = slm(y ~ x + I(x^2) + I(x^3) + sin(2*pi*x) + cos(2*pi*x) + sin(4*pi*x) +
#' cos(4*pi*x) + sin(6*pi*x) + cos(6*pi*x) + sin(8*pi*x) + cos(8*pi*x),
#' method_cov_st = "fitAR", model_selec = -1)
#' summary(reg2)
summary.slm <- function(object, correlation = FALSE, symbolic.cor = FALSE, ...) {
z <- object
p <- z$rank #p = number of variables, with intercept
rdf <- z$df.residual #n-p
#case p==0
if (p == 0) {
r <- z$residuals #hat espilon
n <- length(r) #obs. number
w <- z$weights
if (is.null(w)) {
rss <- sum(r^2)
} else {
rss <- sum(w * r^2)
r <- sqrt(w) * r
}
resvar <- rss/rdf #resvar = sigma2 = rss/n-p = variance
ans <- z[c("call", "terms", if(!is.null(z$weights)) "weights")]
class(ans) <- "summary.slm"
ans$aliased <- is.na(coef(object)) # used in print method
ans$residuals <- r
ans$df <- c(0L, n, length(ans$aliased))
ans$coefficients <- matrix(NA, 0L, 4L)
dimnames(ans$coefficients) <-
list(NULL, c("Estimate", "Std. Error", "z value", "Pr(>|z|)"))
ans$sigma <- sqrt(resvar) #Residual Standard Error
ans$r.squared <- 0
# ans$r.squared <- ans$adj.r.squared <- 0
return(ans)
}
if (is.null(z$terms))
stop("invalid 'slm' object: no 'terms' component")
if(!inherits(object, "slm"))
warning("calling summary.slm(<fake-slm-object>) ...")
Qr <- qr.slm(object)
n <- NROW(Qr$qr) #obs. number
if(is.na(z$df.residual) || n - p != z$df.residual)
warning("residual degrees of freedom in object suggest this is not an \"slm\" fit")
## do not want missing values substituted here
r <- z$residuals #residuals hat epsilon
f <- z$fitted.values #hat Y
w <- z$weights
if (is.null(w)) {
mss <- if (attr(z$terms, "intercept"))
sum((f - mean(f))^2) else sum(f^2)
rss <- sum(r^2)
} else {
mss <- if (attr(z$terms, "intercept")) {
m <- sum(w * f /sum(w))
sum(w * (f - m)^2)
} else sum(w * f^2)
rss <- sum(w * r^2)
r <- sqrt(w) * r
}
resvar <- rss/rdf #sigma2
## see thread at https://stat.ethz.ch/pipermail/r-help/2014-March/367585.html
if (is.finite(resvar) &&
resvar < (mean(f)^2 + var(f)) * 1e-30) # a few times .Machine$double.eps^2
warning("essentially perfect fit: summary may be unreliable")
p1 <- 1L:p #p = var. number
#change table values
norm_matrix = object@norm_matrix
num_var = dim(norm_matrix)[2]
#Cn: covariance matrix estimator
Cn = cov_matrix_estimator(object)
R <- chol2inv(Qr$qr[p1, p1, drop = FALSE]) #= (XtX)^-1
#Modification of the estimation of the standard error of each beta_chapeau_j
se <- sqrt(diag(Cn))/diag(norm_matrix) #sqrt(diag(R) * resvar) #Std. Error of hat beta j
est <- z$coefficients[Qr$pivot[p1]] #Estimate
#Modification of the t-value (-> z-value). Student's test H0 : beta_j = 0.
zval <- est/se #tval <- est/se #"new Student" tests
ans <- z[c("call", "terms", if(!is.null(z$weights)) "weights")]
ans$residuals <- r
ans$coefficients <-
#cbind(est, se, zval, 2*pnorm(abs(zval), mean=0, sd=1, lower.tail=FALSE))
#Modification of the p-value, stars are automatically modified
cbind(est, se, zval, 2*(1-pnorm(abs(zval), mean=0, sd=1))) #2*pt(abs(tval), rdf, lower.tail = FALSE)) #p-value
#attention aux p-values, apparemment ici il faut ecrire 1-pnorm... alors que pour la fstatistic, juste pchisq... et R fait
#le reste automatiquement...?
dimnames(ans$coefficients) <-
list(names(z$coefficients)[Qr$pivot[p1]],
c("Estimate", "Std. Error", "z value", "Pr(>|z|)"))
ans$aliased <- is.na(coef(object)) # used in print method
ans$sigma <- sqrt(resvar)
ans$df <- c(p, rdf, NCOL(Qr$qr)) #p, n-p, p
if (p != attr(z$terms, "intercept")) {
df.int <- if (attr(z$terms, "intercept")) 1L else 0L
ans$r.squared <- mss/(mss + rss)
#ans$adj.r.squared <- 1 - (1 - ans$r.squared) * ((n - df.int)/rdf)
#F-stat: full model vs intercept
Z = solve(expm::sqrtm(as.matrix(Cn[2:num_var,2:num_var])))%*%(norm_matrix[2:num_var,2:num_var]%*%est[-1])
#Z = solve(chol(Cn[2:num_var,2:num_var]))%*%(norm_matrix[2:num_var,2:num_var]%*%est[-1])
ans$chi2statistic <- c(value = sum(Z^2), numdf = num_var - 1)
#ans$fstatistic <- c(value = (mss/(p - df.int))/resvar,
# numdf = p - df.int, dendf = rdf)
} else ans$r.squared <- 0 #else ans$r.squared <- ans$adj.r.squared <- 0
ans$cov.unscaled <- R #cov.unscaled = (XtX)^-1
dimnames(ans$cov.unscaled) <- dimnames(ans$coefficients)[c(1,1)]
if (correlation) {
ans$correlation <- (R * resvar)/outer(se, se) #outer = produit tensoriel
dimnames(ans$correlation) <- dimnames(ans$cov.unscaled)
ans$symbolic.cor <- symbolic.cor
}
if(!is.null(z$na.action)) ans$na.action <- z$na.action
class(ans) <- "summary.slm"
ans
}
#intercept added in the design X. To compute Cn
# #' @title Print summary for slm
# #'
# #' @param x an object of class "\code{summary.slm}", usually, a result of a call to \code{summary.slm}.
# #' @param digits the number of significant digits to use when printing.
# #' @param symbolic.cor logical. If TRUE, print the correlations in a symbolic form (see \code{\link[stats:symnum]{symnum}}) rather than as numbers.
# #' @param signif.stars logical. If TRUE, ‘significance stars’ are printed for each coefficient.
# #' @param ... further arguments passed to or from other methods.
# #'
# #' @import stats
# #'
# #' @seealso
# #' \code{\link[stats:summary.lm]{summary.lm}}.
# #'
#' @export
print.summary.slm <-
function (x, digits = max(3L, getOption("digits") - 3L),
symbolic.cor = x$symbolic.cor,
signif.stars = getOption("show.signif.stars"), ...)
{
cat("\nCall:\n", # S has ' ' instead of '\n'
paste(deparse(x$call), sep="\n", collapse = "\n"), "\n\n", sep = "") #call display
resid <- x$residuals #hat espilon
df <- x$df #n-p
rdf <- df[2L]
cat(if(!is.null(x$weights) && diff(range(x$weights))) "Weighted ",
"Residuals:\n", sep = "")
if (rdf > 5L) {
nam <- c("Min", "1Q", "Median", "3Q", "Max")
rq <- if (length(dim(resid)) == 2L)
structure(apply(t(resid), 1L, quantile),
dimnames = list(nam, dimnames(resid)[[2L]]))
else {
zz <- zapsmall(quantile(resid), digits + 1L)
structure(zz, names = nam)
}
print(rq, digits = digits, ...)
}
else if (rdf > 0L) {
print(resid, digits = digits, ...)
} else { # rdf == 0 : perfect fit!
cat("ALL", df[1L], "residuals are 0: no residual degrees of freedom!")
cat("\n")
}
if (length(x$aliased) == 0L) {
cat("\nNo Coefficients\n")
} else {
if (nsingular <- df[3L] - df[1L])
cat("\nCoefficients: (", nsingular,
" not defined because of singularities)\n", sep = "")
else cat("\nCoefficients:\n")
coefs <- x$coefficients
if(!is.null(aliased <- x$aliased) && any(aliased)) {
cn <- names(aliased)
coefs <- matrix(NA, length(aliased), 4, dimnames=list(cn, colnames(coefs)))
coefs[!aliased, ] <- x$coefficients
}
printCoefmat(coefs, digits = digits, signif.stars = signif.stars,
na.print = "NA", ...)
}
##
cat("\nResidual standard error:",
format(signif(x$sigma, digits)))
# format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom")
if(nzchar(mess <- naprint(x$na.action))) cat(" (",mess, ")\n", sep = "")
if (!is.null(x$chi2statistic)) {
cat("\nMultiple R-squared: ", formatC(x$r.squared, digits = digits))
cat("\nchi2-statistic:", formatC(x$chi2statistic[1L], digits = digits), "on", x$chi2statistic[2L],"DF, p-value:",
format.pval(pchisq(x$chi2statistic[1L], x$chi2statistic[2L], lower.tail = FALSE), digits = digits))
#cat("Multiple R-squared: ", formatC(x$r.squared, digits = digits))
#cat(",\tAdjusted R-squared: ",formatC(x$adj.r.squared, digits = digits),
# "\nF-statistic:", formatC(x$fstatistic[1L], digits = digits),
# "on", x$fstatistic[2L], "and",
# x$fstatistic[3L], "DF, p-value:",
# format.pval(pf(x$fstatistic[1L], x$fstatistic[2L],
# x$fstatistic[3L], lower.tail = FALSE),
# digits = digits))
#cat("\tp-value:", formatC(x$fstatistic[1L], digits = digits))
cat("\n")
}
correl <- x$correlation
if (!is.null(correl)) {
p <- NCOL(correl)
if (p > 1L) {
cat("\nCorrelation of Coefficients:\n")
if(is.logical(symbolic.cor) && symbolic.cor) {# NULL < 1.7.0 objects
print(symnum(correl, abbr.colnames = NULL))
} else {
correl <- format(round(correl, 2), nsmall = 2, digits = digits)
correl[!lower.tri(correl)] <- ""
print(correl[-1, -p, drop=FALSE], quote = FALSE)
}
}
}
cat("\n")#- not in S
invisible(x)
}
# using qr(<lm>) as interface to <lm>$qr :
# #' @title qr.slm
# #'
# #' @param x a numeric or complex matrix whose \code{QR} decomposition is to be computed. Logical matrices are coerced to numeric.
# #' @param ... further arguments passed to or from other methods.
# #'
# #' @seealso
# #' \code{\link[base:qr]{qr}}.
# #'
#' @export
qr.slm <- function(x, ...) {
if(is.null(r <- x$qr))
stop("slm object does not have a proper 'qr' component.
Rank zero or should not have used lm(.., qr=FALSE).")
r
}
#' @title Predict for slm object
#'
#' @description Predicted values based on \code{slm} object.
#'
#' @param object an object of class \code{slm}.
#' @param newdata an optional data frame in which to look for variables with which to predict.
#' \code{newdata} must contain only variables and not the intercept.
#' If omitted, the fitted values are used.
#' @param interval type of interval calculation. It can be only \code{interval = "confidence"}, the default value. It computes
#' the confidence intervals for \eqn{x' beta}, where \eqn{x'} is a new observation of the design.
#' @param level a number between 0 and 1, which indicates the tolerance/confidence level.
#' @param ... further arguments passed to or from other methods.
#'
#' @details This function produces predicted values, obtained by evaluating the regression function in the frame \code{newdata}
#' (which defaults to \code{model.frame(object)}). If \code{newdata} is omitted the predictions are based on the data used for the fit.
#'
#' @return This function produces a vector of predictions or a matrix of predictions and bounds with column names \code{fit}, \code{lwr},
#' and \code{upr} if \code{interval} is set.
#'
#' @importFrom stats predict.lm qnorm lm
#'
#' @seealso
#' \code{\link[stats:predict.lm]{predict.lm}}.
#'
#' @export
#'
#' @examples
#' data("shan")
#' reg1 = slm(shan$PM_Xuhui ~ . , data = shan, method_cov_st = "fitAR", model_selec = -1)
#' predict(reg1)
#'
#' data("co2")
#' y = as.vector(co2)
#' x = as.vector(time(co2)) - 1958
#' reg2 = slm(y ~ x + I(x^2) + I(x^3) + sin(2*pi*x) + cos(2*pi*x) + sin(4*pi*x) +
#' cos(4*pi*x) + sin(6*pi*x) + cos(6*pi*x) + sin(8*pi*x) + cos(8*pi*x),
#' method_cov_st = "fitAR", model_selec = -1)
#' predict(reg2)
predict.slm <- function(object, newdata=NULL, interval="confidence", level=0.95, ...) {
alpha = 1 - level
Dn = object@norm_matrix
Cn = cov_matrix_estimator(object)
design = model.matrix(object)
#design = cbind(rep(1,length(as.matrix(object$model[1]))),as.matrix(object$model[-1]))
tdesign = t(design)
if (interval=="confidence") {
if (is.null(newdata)) {
pred = predict.lm(object, interval = interval, level = level, ...)
#Bhat~N(B,C), x’B~N(xB,x’Cx) IC pour x’B, IC pour une nouvelle observation).
for (i in seq(1,length(pred[,1]))) {
pred[i,2] = pred[i,1] - qnorm(1-(alpha/2))*sqrt(design[i,]%*%solve(Dn)%*%Cn%*%solve(Dn)%*%tdesign[,i])
pred[i,3] = pred[i,1] + qnorm(1-(alpha/2))*sqrt(design[i,]%*%solve(Dn)%*%Cn%*%solve(Dn)%*%tdesign[,i])
}
print(pred)
} else {
mod = object$model
names(mod)[1] = "Y"
reglm = lm(Y~., data=mod)
pred = predict.lm(reglm, newdata, interval = interval, level = level, ...)
newdata = cbind(1,as.matrix(newdata))
tnewdata = t(newdata)
for (i in seq(1,length(pred[,1]))) {
pred[i,2] = pred[i,1] - qnorm(1-(alpha/2))*sqrt(newdata[i,]%*%solve(Dn)%*%Cn%*%solve(Dn)%*%tnewdata[,i])
pred[i,3] = pred[i,1] + qnorm(1-(alpha/2))*sqrt(newdata[i,]%*%solve(Dn)%*%Cn%*%solve(Dn)%*%tnewdata[,i])
}
print(pred)
}
} else if (interval=="prediction") {
print("None available")
} else {
predict.lm(object, newdata = newdata, ...) #none renvoie la prediction
}
}
#' @title Plot.slm
#'
#' @description Same function as the \code{\link[stats:plot.lm]{plot.lm}} function.
#'
#' @param x \code{slm} object.
#' @param ... other parameters to be passed through to plotting functions.
#'
#' @return This function returns the graphics of \code{plot.lm(x)}.
#'
#' @importFrom stats lm
#' @importFrom graphics plot
#'
#' @export
#'
#' @examples
#' data("shan")
#' reg = slm(shan$PM_Xuhui ~ . , data = shan, method_cov_st = "fitAR", model_selec = -1)
#' plot(reg)
plot.slm <- function(x, ...) {
mod = x$model
names(mod)[1] = "Y"
reglm = lm(Y~., data=mod)
plot(reglm, ...)
}
# plot.slm <- function(x, which = c(1:3, 5),
# caption = list("Residuals vs Fitted", "Normal Q-Q",
# "Scale-Location", "Cook's distance",
# "Residuals vs Leverage",
# expression("Cook's dist vs Leverage " * h[ii] / (1 - h[ii]))),
# panel = if(add.smooth) panel.smooth else points,
# sub.caption = NULL, main = "",
# ask = prod(par("mfcol")) < length(which) && dev.interactive(),
# ...,
# id.n = 3, labels.id = names(residuals(x)), cex.id = 0.75,
# qqline = TRUE, cook.levels = c(0.5, 1.0),
# add.smooth = getOption("add.smooth"), label.pos = c(4,2),
# cex.caption = 1, cex.oma.main = 1.25) {
# mod = x$model
# names(mod)[1] = "Y"
# reglm = lm(Y~., data=mod)
#
# plot(reglm, which = c(1:3, 5),
# caption = list("Residuals vs Fitted", "Normal Q-Q",
# "Scale-Location", "Cook's distance",
# "Residuals vs Leverage",
# expression("Cook's dist vs Leverage " * h[ii] / (1 - h[ii]))),
# panel = if(add.smooth) panel.smooth else points,
# sub.caption = NULL, main = "",
# ask = prod(par("mfcol")) < length(which) && dev.interactive(),
# ...,
# id.n = 3, labels.id = names(residuals(x)), cex.id = 0.75,
# qqline = TRUE, cook.levels = c(0.5, 1.0),
# add.smooth = getOption("add.smooth"), label.pos = c(4,2),
# cex.caption = 1, cex.oma.main = 1.25)
#
# }
#' Calculate Variance-Covariance Matrix for a Fitted Model Object of class slm
#'
#' @description Returns the variance-covariance matrix of the (non-normalized) least squares estimators for an object of class slm.
#'
#' @param object a fitted model object of class slm.
#' @param ... additional arguments for method functions.
#'
#' @return The variance-covariance matrix of the (non-normalized) least squares estimators for an object of class slm.
#'
#' @importFrom stats vcov
#'
#' @export
#'
#' @seealso
#' The generic function \code{\link[stats:vcov]{vcov}}.
#'
#' The function \code{\link[slm:cov_matrix_estimator]{cov_matrix_estimator}}.
#'
#' @examples
#' n = 500
#' eps = generative_process(n,"AR1",c(0.7))
#' X = as.matrix(generative_model(n,"mod2"))
#' Y = 3 + 2*X[,2] + eps
#' reg = slm(Y ~ X, method_cov_st = "fitAR", model_selec = -1)
#' vcov(reg)
vcov.slm <- function(object, ...) {
vc = vcov(lm(object), ...)
vcslm = solve(object@norm_matrix)%*%cov_matrix_estimator(object)%*%solve(object@norm_matrix)
dimnames(vcslm) = dimnames(vc)
return(vcslm)
#Dn^-1*Cn*Dn^-1, enlever la normalisation de Cn
}
#vcov = matrice de cov des beta chapeaux non normalises (sigma 2 XtX-1 dans le cas iid).
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.