當前位置: 首頁>>編程示例 >>用法及示例精選 >>正文


R shash Sinh-arcsinh 位置比例和形狀模型族

R語言 shash 位於 mgcv 包(package)。

說明

shash 係列實現了 Jones 和 Pewsey (2009) 的 four-parameter sinh-arcsinh(哈希)分布。密度的位置、尺度、偏度和峰度可以取決於加性平滑預測因子。僅可與 gam 一起使用,線性預測變量通過公式列表指定。在使用這種靈活模型之前,值得仔細考慮數據是否足以支持估計。

用法

shash(link = list("identity", "logeb", "identity", "identity"), 
      b = 1e-2, phiPen = 1e-3)

參數

link

由四個字符組成的向量,指示位置、尺度、偏度和峰度參數的鏈接函數。

b

logeb 鏈接函數的正參數,請參閱詳細信息。

phiPen

峰度參數的脊罰分的正乘數。除非您知道自己在做什麽,否則請勿觸摸它,請參閱詳細信息。

細節

shash族的密度函數是

其中 。這裏 分別控製位置和比例, 確定偏度,而 控製尾重。 shash 可以對任一側的偏度進行建模,具體取決於 的符號。此外,shash 的尾部可以比普通尾部更輕 ( ) 或更重 ( )。出於擬合目的,這裏我們使用

密度基於第 4.1 節第二行給出的表達式以及 Jones 和 Pewsey (2009) 的方程 (2),並使用第 4.3 節給出的簡單重新參數化。

用於 的鏈接函數是 logeb ,其中是 ,因此反向鏈接是 。關鍵是我們不允許 變得小於一個小常數 b。可能性包括嶺罰分 ,它將 縮小到零。當有足夠的數據可用時,嶺罰分不會對擬合產生太大影響,但在將模型擬合到小數據集時將其包含在內很有用,以避免 發散到無窮大(Jones 和 Pewsey 已經確定了這個問題( 2009))。

繼承自類 general.family 的對象。

例子


###############
# Shash dataset
###############
##  Simulate some data from shash
set.seed(847)
n <- 1000
x <- seq(-4, 4, length.out = n)

X <- cbind(1, x, x^2)
beta <- c(4, 1, 1)
mu <- X %*% beta 

sigma =  .5+0.4*(x+4)*.5            # Scale
eps = 2*sin(x)                      # Skewness
del = 1 + 0.2*cos(3*x)              # Kurtosis

dat <-  mu + (del*sigma)*sinh((1/del)*asinh(qnorm(runif(n))) + (eps/del))
dataf <- data.frame(cbind(dat, x))
names(dataf) <- c("y", "x")
plot(x, dat, xlab = "x", ylab = "y")

## Fit model
fit <- gam(list(y ~ s(x), # <- model for location 
                  ~ s(x),   # <- model for log-scale
                  ~ s(x),   # <- model for skewness
                  ~ s(x, k = 20)), # <- model for log-kurtosis
           data = dataf, 
           family = shash, # <- new family 
           optimizer = "efs")

## Plotting truth and estimates for each parameters of the density 
muE <- fit$fitted[ , 1]
sigE <- exp(fit$fitted[ , 2])
epsE <- fit$fitted[ , 3]
delE <- exp(fit$fitted[ , 4])

par(mfrow = c(2, 2))
plot(x, muE, type = 'l', ylab = expression(mu(x)), lwd = 2)
lines(x, mu, col = 2, lty = 2, lwd = 2)
legend("top", c("estimated", "truth"), col = 1:2, lty = 1:2, lwd = 2)

plot(x, sigE, type = 'l', ylab = expression(sigma(x)), lwd = 2)
lines(x, sigma, col = 2, lty = 2, lwd = 2)

plot(x, epsE, type = 'l', ylab = expression(epsilon(x)), lwd = 2)
lines(x, eps, col = 2, lty = 2, lwd = 2)

plot(x, delE, type = 'l', ylab = expression(delta(x)), lwd = 2)
lines(x, del, col = 2, lty = 2, lwd = 2)

## Plotting true and estimated conditional density
par(mfrow = c(1, 1))
plot(x, dat, pch = '.', col = "grey", ylab = "y", ylim = c(-35, 70))
for(qq in c(0.001, 0.01, 0.1, 0.5, 0.9, 0.99, 0.999)){
  est <- fit$family$qf(p=qq, mu = fit$fitted)
  true <- mu + (del * sigma) * sinh((1/del) * asinh(qnorm(qq)) + (eps/del))
  lines(x, est, type = 'l', col = 1, lwd = 2)
  lines(x, true, type = 'l', col = 2, lwd = 2, lty = 2)
}
legend("topleft", c("estimated", "truth"), col = 1:2, lty = 1:2, lwd = 2)

#####################
## Motorcycle example
#####################

# Here shash is overkill, in fact the fit is not good, relative
# to what we would get with mgcv::gaulss
library(MASS)

b <- gam(list(accel~s(times, k=20, bs = "ad"), ~s(times, k = 10), ~1, ~1),
         data=mcycle, family=shash)

par(mfrow = c(1, 1))
xSeq <- data.frame(cbind("accel" = rep(0, 1e3),
                   "times" = seq(2, 58, length.out = 1e3)))
pred <- predict(b, newdata = xSeq)
plot(mcycle$times, mcycle$accel, ylim = c(-180, 100))
for(qq in c(0.1, 0.3, 0.5, 0.7, 0.9)){
  est <- b$family$qf(p=qq, mu = pred)
  lines(xSeq$times, est, type = 'l', col = 2)
}

plot(b, pages = 1, scale = FALSE)

作者

Matteo Fasiolo <matteo.fasiolo@gmail.com> and Simon N. Wood.

參考

Jones, M. and A. Pewsey (2009). Sinh-arcsinh distributions. Biometrika 96 (4), 761-780. doi:10.1093/biomet/asp053

Wood, S.N., N. Pya and B. Saefken (2016), Smoothing parameter and model selection for general smooth models. Journal of the American Statistical Association 111, 1548-1575 doi:10.1080/01621459.2016.1180986

相關用法


注:本文由純淨天空篩選整理自R-devel大神的英文原創作品 Sinh-arcsinh location scale and shape model family。非經特殊聲明,原始代碼版權歸原作者所有,本譯文未經允許或授權,請勿轉載或複製。