我正在使用一个小数据集(21 个观察值),并且在 R 中有以下正常的 QQ 图:
看到该图不支持正态性,我可以推断出基本分布吗?在我看来,更偏向右侧的分布会更合适,对吗?此外,我们还能从数据中得出什么其他结论?
我正在使用一个小数据集(21 个观察值),并且在 R 中有以下正常的 QQ 图:
看到该图不支持正态性,我可以推断出基本分布吗?在我看来,更偏向右侧的分布会更合适,对吗?此外,我们还能从数据中得出什么其他结论?
如果这些值沿着一条线分布,则分布具有与我们假设的理论分布相同的形状(直至位置和比例)。
局部行为:当查看 y 轴上的排序样本值和 x 轴上的(近似)预期分位数时,我们可以通过查看是否值或多或少比理论分布在该图的该部分中假设的集中:
正如我们所看到的,较不集中的点增加越来越多,集中点的增加速度比整体线性关系所暗示的要慢,并且在极端情况下对应于样本密度的差距(显示为接近垂直的跳跃)或常量值的峰值(水平对齐的值)。这使我们能够发现重尾或轻尾,因此偏度大于或小于理论分布,依此类推。
整体外观:
这是 QQ 图的平均样子(对于特定的分布选择):
但是随机性往往会使事情变得模糊,尤其是对于小样本:
请注意,在时,结果可能比那里显示的变化更大——我生成了几组这样的六个图,并选择了一个“不错”的集合,您可以在其中同时看到所有六个图中的形状。有时直的关系看起来是弯曲的,弯曲的关系看起来是直的,重尾看起来是歪斜的,等等——对于这么小的样本,情况通常可能不太清楚:
可以识别出比那些特征更多的特征(例如离散性),但是在的情况下,即使是这样的基本特征也可能难以发现;我们不应该试图“过度解释”每一个小小的摆动。随着样本量变大,一般来说,图“稳定”并且特征变得更清晰可解释而不是表示噪声。[对于一些非常重的分布,即使在相当大的样本量下,罕见的大异常值也可能会阻止图片很好地稳定。]
当您尝试决定您应该担心多少曲率或摆动时,您可能还会发现此处的建议很有用。
一般而言,更适合解释的指南还包括以更小和更大的样本大小显示。
我做了一个闪亮的应用程序来帮助解释正常的 QQ 情节。试试这个链接。
在这个应用程序中,您可以调整数据的偏度、尾度(峰度)和模态,您可以看到直方图和 QQ 图的变化情况。相反,您可以按照给定 QQ 图的模式来使用它,然后检查偏度等应该如何。
有关详细信息,请参阅其中的文档。
我意识到我没有足够的可用空间来在线提供这个应用程序。根据要求,我将提供所有三个代码块:sample.R
和server.R
这里ui.R
。有兴趣运行此应用程序的人可以将这些文件加载到 Rstudio 中,然后在您自己的 PC 上运行它。
sample.R
文件:
# Compute the positive part of a real number x,
# which is $\max(x, 0)$.
positive_part <- function(x) {ifelse(x > 0, x, 0)}
# This function generates n data points from some
# unimodal population.
# Input: ----------------------------------------------------
# n: sample size;
# mu: the mode of the population, default value is 0.
# skewness: the parameter that reflects the skewness of the
# distribution, note it is not
# the exact skewness defined in statistics textbook,
# the default value is 0.
# tailedness: the parameter that reflects the tailedness
# of the distribution, note it is
# not the exact kurtosis defined in textbook,
# the default value is 0.
# When all arguments take their default values, the data will
# be generated from standard
# normal distribution.
random_sample <- function(n, mu = 0, skewness = 0, §
tailedness = 0){
sigma = 1
# The sampling scheme resembles the rejection sampling.
# For each step, an initial data point
# was proposed, and it will be rejected or accepted based on
# the weights determined by the
# skewness and tailedness of input.
reject_skewness <- function(x){
scale = 1
# if `skewness` > 0 (means data are right-skewed),
# then small values of x will be rejected
# with higher probability.
l <- exp(-scale * skewness * x)
l/(1 + l)
}
reject_tailedness <- function(x){
scale = 1
# if `tailedness` < 0 (means data are lightly-tailed),
# then big values of x will be rejected with
# higher probability.
l <- exp(-scale * tailedness * abs(x))
l/(1 + l)
}
# w is another layer option to control the tailedness, the
# higher the w is, the data will be
# more heavily-tailed.
w = positive_part((1 - exp(-0.5 * tailedness)))/(1 +
exp(-0.5 * tailedness))
filter <- function(x){
# The proposed data points will be accepted only if it
# satified the following condition,
# in which way we controlled the skewness and tailedness of
# data. (For example, the
# proposed data point will be rejected more frequently if it
# has higher skewness or
# tailedness.)
accept <- runif(length(x)) > reject_tailedness(x) *
reject_skewness(x)
x[accept]
}
result <- filter(mu + sigma * ((1 - w) * rnorm(n) + w * rt(n, 5)))
# Keep generating data points until the length of data vector
# reaches n.
while (length(result) < n) {
result <- c(result, filter(mu + sigma * ((1 - w) * rnorm(n) +
w * rt(n, 5))))
}
result[1:n]
}
multimodal <- function(n, Mu, skewness = 0, tailedness = 0) {
# Deal with the bimodal case.
mumu <- as.numeric(Mu %*% rmultinom(n, 1, rep(1, length(Mu))))
mumu + random_sample(n, skewness = skewness,
tailedness = tailedness)
}
server.R
文件:
library(shiny)
# Need 'ggplot2' package to get a better aesthetic effect.
library(ggplot2)
# The 'sample.R' source code is used to generate data to be
# plotted, based on the input skewness,
# tailedness and modality. For more information, see the source
# code in 'sample.R' code.
source("sample.R")
shinyServer(function(input, output) {
# We generate 10000 data points from the distribution which
# reflects the specification of skewness,
# tailedness and modality.
n = 10000
# 'scale' is a parameter that controls the skewness and
# tailedness.
scale = 1000
# The `reactive` function is a trick to accelerate the app,
# which enables us only generate the data
# once to plot two plots. The generated sample was stored in
# the `data` object to be called later.
data <- reactive({
# For `Unimodal` choice, we fix the mode at 0.
if (input$modality == "Unimodal") {mu = 0}
# For `Bimodal` choice, we fix the two modes at -2 and 2.
if (input$modality == "Bimodal") {mu = c(-2, 2)}
# Details will be explained in `sample.R` file.
sample1 <- multimodal(n, mu, skewness = scale *
input$skewness, tailedness = scale * input$kurtosis)
data.frame(x = sample1)})
output$histogram <- renderPlot({
# Plot the histogram.
ggplot(data(), aes(x = x)) +
geom_histogram(aes(y = ..density..), binwidth = .5,
colour = "black", fill = "white") +
xlim(-6, 6) +
# Overlay the density curve.
geom_density(alpha = .5, fill = "blue") +
ggtitle("Histogram of Data") +
theme(plot.title = element_text(lineheight = .8,
face = "bold"))
})
output$qqplot <- renderPlot({
# Plot the QQ plot.
ggplot(data(), aes(sample = x)) + stat_qq() +
ggtitle("QQplot of Data") +
theme(plot.title = element_text(lineheight=.8,
face = "bold"))
})
})
最后,ui.R
文件:
library(shiny)
# Define UI for application that helps students interpret the
# pattern of (normal) QQ plots.
# By using this app, we can show students the different patterns
# of QQ plots (and the histograms,
# for completeness) for different type of data distributions.
# For example, left skewed heavy tailed
# data, etc.
# This app can be (and is encouraged to be) used in a reversed
# way, namely, show the QQ plot to the
# students first, then tell them based on the pattern of the QQ
# plot, the data is right skewed, bimodal,
# heavy-tailed, etc.
shinyUI(fluidPage(
# Application title
titlePanel("Interpreting Normal QQ Plots"),
sidebarLayout(
sidebarPanel(
# The first slider can control the skewness of input data.
# "-1" indicates the most left-skewed
# case while "1" indicates the most right-skewed case.
sliderInput("skewness", "Skewness", min = -1, max = 1,
value = 0, step = 0.1, ticks = FALSE),
# The second slider can control the skewness of input data.
# "-1" indicates the most light tail
# case while "1" indicates the most heavy tail case.
sliderInput("kurtosis", "Tailedness", min = -1, max = 1,
value = 0, step = 0.1, ticks = FALSE),
# This selectbox allows user to choose the number of modes
# of data, two options are provided:
# "Unimodal" and "Bimodal".
selectInput("modality", label = "Modality",
choices = c("Unimodal" = "Unimodal",
"Bimodal" = "Bimodal"),
selected = "Unimodal"),
br(),
# The following helper information will be shown on the
# user interface to give necessary
# information to help users understand sliders.
helpText(p("The skewness of data is controlled by moving
the", strong("Skewness"), "slider,",
"the left side means left skewed while the right
side means right skewed."),
p("The tailedness of data is controlled by moving
the", strong("Tailedness"), "slider,",
"the left side means light tailed while the
right side means heavy tailed."),
p("The modality of data is controlled by selecting
the modality from", strong("Modality"),
"select box.")
)
),
# The main panel outputs two plots. One plot is the histogram
# of data (with the non-parametric density
# curve overlaid), to get a better visualization, we restricted
# the range of x-axis to -6 to 6 so
# that part of the data will not be shown when heavy-tailed
# input is chosen. The other plot is the
# QQ plot of data, as convention, the x-axis is the theoretical
# quantiles for standard normal distri-
# bution and the y-axis is the sample quantiles of data.
mainPanel(
plotOutput("histogram"),
plotOutput("qqplot")
)
)
)
)
教授给出了一个非常有用(且直观)的解释。Philippe Rigollet 参加 MIT MOOC 课程:18.650 Statistics for Applications,2016 年秋季 - 观看 45 分钟的视频
https://www.youtube.com/watch?v=vMaKx9fmJHE
我粗略地复制了他的图表,我把它放在我的笔记中,因为我发现它非常有用。
在示例 1 中,在左上图中,我们看到右尾的经验(或样本)分位数小于理论分位数
Qe < Qt
这可以使用概率密度函数来解释。对于相同的值,经验分位数在理论分位数的左边,这意味着经验分布的右尾比理论分布的右尾“轻”,即它更快地下降到接近零的值。
由于该线程被认为是权威的“如何解释正常 qq 图”StackExchange 帖子,我想向读者指出正常 qq 图和超峰态统计之间的良好、精确的数学关系。
这里是:
https://stats.stackexchange.com/a/354076/102879
一个简短(而且过于简化)的总结如下(有关更精确的数学陈述,请参见链接):您实际上可以将正常 qq 图中的过度峰度视为数据分位数与相应的理论正态分位数之间的平均距离,加权从数据到平均值的距离。因此,当 qq 图尾部的绝对值通常在极端方向上大大偏离预期的正常值时,您就有正的超峰度。
因为峰度是这些偏差的平均值,由与平均值的距离加权,所以靠近 qq 图中心的值对峰度几乎没有影响。因此,过度峰度与分布的中心无关,即“峰值”所在的位置。相反,过度峰态几乎完全由数据分布尾部与正态分布的比较决定。