In the previous post, we examined the Shapiro-Wilk normality test and the familiar R convention of using four functions for a probability distribution:
| Prefix | Meaning | Normal example |
|---|---|---|
| d | Density | dnorm() |
| p | Cumulative probability | pnorm() |
| q | Quantile | qnorm() |
| r | Random generation | rnorm() |
But Base R does not provide:
dshapiro()
pshapiro()
qshapiro()
rshapiro()
Why not?
The reason is important. The Shapiro-Wilk statistic W does not have a simple probability distribution such as the normal, t, chi-square or F distribution. Its null distribution depends strongly on the sample size n, and the p-values reported by shapiro.test() are obtained using numerical approximations.
That does not stop us from building our own educational d, p, q and r functions.
We can approximate the null distribution of W using Monte Carlo simulation.
1. What distribution are we actually studying?
Suppose:
For each random sample, calculate the Shapiro-Wilk statistic:
If we repeat this experiment thousands of times, we obtain a distribution of W under the null hypothesis.
That distribution is:
The sampling distribution of the Shapiro-Wilk W statistic when the data are truly normal.
This is the distribution we will use to construct:
dshapiro()= approximate density of W;pshapiro()= approximate cumulative probability for W;qshapiro()= approximate quantiles of W;rshapiro()= random values of W under normality.
2. A crucial fact: the distribution depends on sample size
Unlike the standard normal distribution, the distribution of W changes with n.
For example, the distribution of W for:
n = 30
n = 100
will be different.
Therefore every function we write must include n as an argument.
Conceptually:
pshapiro(w, n = 30)
means:
For samples of size 30 drawn from a normal population, what proportion of Shapiro-Wilk W statistics are less than or equal to w?
3. First build rshapiro()
The most natural function to build first is rshapiro().
To generate one random Shapiro-Wilk W value:
- generate n observations from a normal distribution;
- run
shapiro.test(); - extract W.
Then repeat this procedure as many times as requested.
4. Base R implementation of rshapiro()
rshapiro <- function(
N,
n,
mean = 0,
sd = 1
) {
# N = number of W statistics to generate
# n = sample size used for each Shapiro-Wilk test
if (N < 1)
return(numeric(0))
if (n < 3)
stop("Shapiro-Wilk requires at least 3 observations")
if (n > 5000)
stop("Base R shapiro.test() allows at most 5000 observations")
if (sd <= 0)
stop("sd must be positive")
# Create storage for simulated W statistics
W <- numeric(N)
# Repeat N times
for (i in seq_len(N)) {
# Generate a sample from a normal population
x <- rnorm(
n = n,
mean = mean,
sd = sd
)
# Calculate Shapiro-Wilk statistic
W[i] <- unname(
shapiro.test(x)$statistic
)
}
return(W)
}
5. Demonstrating rshapiro()
Generate ten Shapiro-Wilk statistics for samples of size 20:
set.seed(123)
rshapiro(
N = 10,
n = 20
)
You might obtain values such as:
0.963
0.974
0.952
0.983
0.968
...
Notice that most values are close to 1 because the underlying samples were generated from a normal distribution.
6. Generate the complete null distribution
For teaching purposes, generate a much larger number of W statistics:
set.seed(123)
W30 <- rshapiro(
N = 10000,
n = 30
)
Now inspect the distribution:
hist(
W30,
breaks = 50,
probability = TRUE,
main = "Null distribution of Shapiro-Wilk W",
xlab = "W"
)
This is an approximation to the sampling distribution of W for:
7. Why the mean and standard deviation do not matter
Try:
set.seed(123)
W1 <- rshapiro(
N = 5000,
n = 30,
mean = 0,
sd = 1
)
set.seed(123)
W2 <- rshapiro(
N = 5000,
n = 30,
mean = 100,
sd = 25
)
Compare:
mean(W1)
mean(W2)
sd(W1)
sd(W2)
They will be essentially identical.
This illustrates an important property of the Shapiro-Wilk test:
The W statistic is invariant to changes in location and scale.
Therefore:
is sufficient when simulating the null distribution.
8. Building pshapiro()
The p-function for a distribution gives:
For example:
pshapiro(
q = 0.95,
n = 30
)
asks:
Under normality, what fraction of samples of size 30 produce W less than or equal to 0.95?
9. Monte Carlo implementation of pshapiro()
pshapiro <- function(
q,
n,
B = 10000,
lower.tail = TRUE,
seed = NULL
) {
# q = observed W value
# n = sample size
# B = number of Monte Carlo simulations
if (!is.null(seed))
set.seed(seed)
if (n < 3)
stop("n must be at least 3")
if (n > 5000)
stop("n cannot exceed 5000")
# Generate null distribution of W
Wnull <- rshapiro(
N = B,
n = n
)
# Empirical cumulative probability
p <- sapply(
q,
function(value) {
mean(Wnull <= value)
}
)
# Return lower or upper tail
if (!lower.tail)
p <- 1 - p
return(p)
}
10. Demonstrating pshapiro()
set.seed(100)
pshapiro(
q = 0.95,
n = 30,
B = 20000
)
Suppose this returns approximately:
0.17
This means:
About 17% of normally distributed samples of size 30 produce a Shapiro-Wilk W statistic of 0.95 or smaller.
11. Why Shapiro-Wilk uses the lower tail
For Shapiro-Wilk:
Therefore the p-value is based on:
rather than the upper tail.
This is different from many familiar test statistics such as chi-square or F, where large values provide evidence against the null hypothesis.
12. Compare pshapiro() with shapiro.test()
Generate a sample:
set.seed(42)
x <- rnorm(30)
test <- shapiro.test(x)
test
Extract W:
Wobs <- unname(
test$statistic
)
Wobs
Now estimate its p-value by simulation:
pshapiro(
q = Wobs,
n = length(x),
B = 50000,
seed = 100
)
Compare with:
test$p.value
The values should be reasonably similar, although they will not be identical.
Why?
Our function uses Monte Carlo simulation, whereas shapiro.test() uses specialised analytical approximations.
13. A better Monte Carlo p-value
A simple simulation can occasionally return p = 0 if no simulated W is as small as the observed value.
A common finite-simulation correction is:
We can modify pshapiro() accordingly.
pshapiro <- function(
q,
n,
B = 10000,
lower.tail = TRUE,
seed = NULL
) {
if (!is.null(seed))
set.seed(seed)
if (n < 3)
stop("n must be at least 3")
if (n > 5000)
stop("n cannot exceed 5000")
Wnull <- rshapiro(
N = B,
n = n
)
p <- sapply(
q,
function(value) {
lower <-
(sum(Wnull <= value) + 1) /
(B + 1)
if (lower.tail)
return(lower)
return(1 - lower)
}
)
return(p)
}
14. Building qshapiro()
The q-function is the inverse of the p-function.
It answers:
What value of W corresponds to a particular cumulative probability?
For example:
qshapiro(
p = 0.05,
n = 30
)
asks:
What is approximately the lower 5th percentile of W when n = 30 and the data are normal?
This value acts as an approximate 5% critical value.
15. Base R implementation of qshapiro()
qshapiro <- function(
p,
n,
B = 10000,
lower.tail = TRUE,
seed = NULL
) {
if (!is.null(seed))
set.seed(seed)
if (any(p < 0 | p > 1))
stop("p must lie between 0 and 1")
if (n < 3)
stop("n must be at least 3")
if (n > 5000)
stop("n cannot exceed 5000")
# Generate null distribution
Wnull <- rshapiro(
N = B,
n = n
)
# Convert upper-tail probabilities if requested
probs <- p
if (!lower.tail)
probs <- 1 - p
# Empirical quantiles
result <- quantile(
Wnull,
probs = probs,
names = FALSE,
type = 8
)
return(result)
}
16. Demonstrating qshapiro()
set.seed(123)
qshapiro(
p = 0.05,
n = 30,
B = 50000
)
Suppose this gives:
0.93
Then approximately 5% of normal samples of size 30 generate:
Thus 0.93 is an approximate 5% critical value.
17. Demonstrating the p-q relationship
Suppose:
crit <- qshapiro(
p = 0.05,
n = 30,
B = 50000,
seed = 10
)
crit
Now evaluate its cumulative probability:
pshapiro(
q = crit,
n = 30,
B = 50000,
seed = 20
)
The answer should be close to:
0.05
This illustrates the same relationship we see with:
pnorm(qnorm(0.05))
18. Building dshapiro()
The d-function is slightly more complicated.
For a continuous distribution:
But we do not have a closed-form density equation for W.
Instead, we can:
- simulate thousands of W values;
- estimate their density using kernel density estimation;
- interpolate that density at the requested W value.
19. Base R implementation of dshapiro()
dshapiro <- function(
x,
n,
B = 20000,
bw = "nrd0",
seed = NULL
) {
if (!is.null(seed))
set.seed(seed)
if (n < 3)
stop("n must be at least 3")
if (n > 5000)
stop("n cannot exceed 5000")
# Simulate null distribution
Wnull <- rshapiro(
N = B,
n = n
)
# Kernel density estimate
dens <- density(
Wnull,
bw = bw,
from = 0,
to = 1
)
# Interpolate density at requested x values
answer <- approx(
x = dens$x,
y = dens$y,
xout = x,
rule = 2
)$y
# W cannot theoretically be outside 0 to 1
answer[x < 0 | x > 1] <- 0
return(answer)
}
20. Demonstrating dshapiro()
dshapiro(
x = 0.95,
n = 30,
B = 50000,
seed = 123
)
This gives an estimated probability density around:
Remember:
A density is not a probability.
For a continuous statistic:
The density describes how concentrated the probability distribution is near 0.95.
21. Plotting our estimated Shapiro-Wilk density
wgrid <- seq(
0.80,
1,
length.out = 200
)
dens30 <- dshapiro(
x = wgrid,
n = 30,
B = 50000,
seed = 123
)
plot(
wgrid,
dens30,
type = "l",
lwd = 2,
xlab = "Shapiro-Wilk W",
ylab = "Estimated density",
main = "Estimated null density of W, n = 30"
)
22. Comparing different sample sizes
This is one of the most informative demonstrations.
wgrid <- seq(
0.75,
1,
length.out = 300
)
d10 <- dshapiro(
wgrid,
n = 10,
B = 30000,
seed = 1
)
d30 <- dshapiro(
wgrid,
n = 30,
B = 30000,
seed = 2
)
d100 <- dshapiro(
wgrid,
n = 100,
B = 30000,
seed = 3
)
plot(
wgrid,
d10,
type = "l",
lwd = 2,
xlab = "W",
ylab = "Estimated density",
main = "Null distribution of W depends on sample size"
)
lines(
wgrid,
d30,
lwd = 2,
lty = 2
)
lines(
wgrid,
d100,
lwd = 2,
lty = 3
)
legend(
"topleft",
legend = c(
"n = 10",
"n = 30",
"n = 100"
),
lty = c(1, 2, 3),
lwd = 2
)
As n increases, the distribution of W becomes increasingly concentrated near 1.
This is why a value such as:
may be completely ordinary for a small sample but unusually low for a very large sample.
23. Collecting all four functions together
Here are the complete educational functions in one block.
############################################################
# rshapiro()
############################################################
rshapiro <- function(
N,
n
) {
if (N < 1)
return(numeric(0))
if (n < 3)
stop("n must be at least 3")
if (n > 5000)
stop("n cannot exceed 5000")
W <- numeric(N)
for (i in seq_len(N)) {
x <- rnorm(n)
W[i] <-
unname(
shapiro.test(x)$statistic
)
}
return(W)
}
############################################################
# pshapiro()
############################################################
pshapiro <- function(
q,
n,
B = 10000,
lower.tail = TRUE,
seed = NULL
) {
if (!is.null(seed))
set.seed(seed)
Wnull <- rshapiro(
N = B,
n = n
)
p <- sapply(
q,
function(value) {
lower <-
(sum(Wnull <= value) + 1) /
(B + 1)
if (lower.tail)
return(lower)
return(1 - lower)
}
)
return(p)
}
############################################################
# qshapiro()
############################################################
qshapiro <- function(
p,
n,
B = 10000,
lower.tail = TRUE,
seed = NULL
) {
if (!is.null(seed))
set.seed(seed)
if (any(p < 0 | p > 1))
stop("p must lie between 0 and 1")
Wnull <- rshapiro(
N = B,
n = n
)
probs <- p
if (!lower.tail)
probs <- 1 - p
quantile(
Wnull,
probs = probs,
names = FALSE,
type = 8
)
}
############################################################
# dshapiro()
############################################################
dshapiro <- function(
x,
n,
B = 20000,
bw = "nrd0",
seed = NULL
) {
if (!is.null(seed))
set.seed(seed)
Wnull <- rshapiro(
N = B,
n = n
)
dens <- density(
Wnull,
bw = bw,
from = 0,
to = 1
)
answer <- approx(
x = dens$x,
y = dens$y,
xout = x,
rule = 2
)$y
answer[
x < 0 |
x > 1
] <- 0
return(answer)
}
24. Demonstrating all four functions
############################################################
# Choose sample size
############################################################
n <- 30
############################################################
# r: generate random W values
############################################################
set.seed(1)
rshapiro(
N = 5,
n = n
)
############################################################
# p: cumulative probability
############################################################
pshapiro(
q = 0.95,
n = n,
B = 30000,
seed = 1
)
############################################################
# q: lower 5% critical value
############################################################
qshapiro(
p = 0.05,
n = n,
B = 30000,
seed = 1
)
############################################################
# d: estimated density
############################################################
dshapiro(
x = 0.95,
n = n,
B = 30000,
seed = 1
)
25. Reproducing a complete Shapiro-Wilk test manually
Now we can use our functions to construct an approximate Shapiro-Wilk test.
set.seed(123)
x <- rnorm(
30,
mean = 50,
sd = 10
)
# Calculate observed W
Wobs <-
unname(
shapiro.test(x)$statistic
)
Wobs
Now obtain our Monte Carlo p-value:
p_manual <- pshapiro(
q = Wobs,
n = length(x),
B = 50000,
seed = 999
)
p_manual
Compare with R:
shapiro.test(x)$p.value
The two values should usually be close.
26. Writing our own complete my_shapiro_test()
my_shapiro_test <- function(
x,
B = 50000,
seed = NULL
) {
# Remove missing values
x <- x[!is.na(x)]
n <- length(x)
if (n < 3)
stop("At least 3 observations are required")
if (n > 5000)
stop("At most 5000 observations are allowed")
# Calculate observed Shapiro-Wilk statistic
Wobs <-
unname(
shapiro.test(x)$statistic
)
# Estimate lower-tail probability
pvalue <- pshapiro(
q = Wobs,
n = n,
B = B,
seed = seed
)
result <- list(
statistic = Wobs,
p.value = pvalue,
sample.size = n,
simulations = B
)
return(result)
}
Example:
set.seed(42)
x <- rnorm(30)
my_shapiro_test(
x,
B = 50000,
seed = 100
)
Compare:
shapiro.test(x)
27. Test our implementation on non-normal data
Generate exponential observations:
set.seed(50)
x <- rexp(
30,
rate = 1
)
my_shapiro_test(
x,
B = 50000,
seed = 100
)
shapiro.test(x)
Both procedures should generally indicate strong evidence against normality.
28. Visualising the observed statistic inside the null distribution
This is perhaps the best way to understand a hypothesis test.
set.seed(50)
x <- rexp(30)
Wobs <-
unname(
shapiro.test(x)$statistic
)
set.seed(100)
Wnull <- rshapiro(
N = 50000,
n = length(x)
)
hist(
Wnull,
breaks = 60,
probability = TRUE,
main = "Observed W within the null distribution",
xlab = "Shapiro-Wilk W"
)
abline(
v = Wobs,
lwd = 3
)
For strongly non-normal data, the observed W will sit deep in the lower tail.
The p-value is simply the fraction of the null distribution lying at or below that observed value.
29. Mark the rejection region
critical <-
quantile(
Wnull,
0.05
)
critical
Then:
hist(
Wnull,
breaks = 60,
main = "Shapiro-Wilk null distribution",
xlab = "W"
)
abline(
v = critical,
lwd = 2,
lty = 2
)
abline(
v = Wobs,
lwd = 3
)
The rejection rule is approximately:
30. Understanding d, p, q and r through Shapiro-Wilk
We can now express the four functions conceptually.
rshapiro()
Generates random Shapiro-Wilk statistics:
pshapiro()
Calculates:
qshapiro()
Finds w such that:
dshapiro()
Estimates the local probability density:
31. Comparing with the normal d-p-q-r functions
| Purpose | Normal distribution | Shapiro-Wilk W distribution |
|---|---|---|
| Density | dnorm() |
dshapiro() |
| Cumulative probability | pnorm() |
pshapiro() |
| Quantile | qnorm() |
qshapiro() |
| Random generation | rnorm() |
rshapiro() |
32. One important difference
For the normal distribution:
dnorm()
pnorm()
qnorm()
rnorm()
are based on well-established numerical algorithms for a mathematically defined probability distribution.
Our:
dshapiro()
pshapiro()
qshapiro()
rshapiro()
are Monte Carlo approximations.
Their precision therefore depends on:
- the number of simulations B;
- the random seed;
- the kernel bandwidth used for density estimation.
33. Monte Carlo error
Suppose the true probability is p.
If we estimate it from B simulations, the Monte Carlo standard error is approximately:
For:
B = 10000
the standard error is approximately:
Increasing B reduces simulation noise.
34. Demonstrating Monte Carlo variability
pshapiro(
0.95,
n = 30,
B = 1000,
seed = 1
)
pshapiro(
0.95,
n = 30,
B = 1000,
seed = 2
)
pshapiro(
0.95,
n = 30,
B = 1000,
seed = 3
)
The answers will vary slightly.
Now increase B:
pshapiro(
0.95,
n = 30,
B = 100000,
seed = 1
)
The estimate becomes more stable.
35. A faster approach: simulate once and reuse
The functions above regenerate the null distribution every time. This is useful pedagogically, but inefficient.
For repeated calculations, generate the distribution once:
set.seed(123)
Wnull30 <- rshapiro(
N = 100000,
n = 30
)
Then use it directly.
p-function
mean(
Wnull30 <= 0.95
)
q-function
quantile(
Wnull30,
0.05
)
d-function
d <- density(
Wnull30,
from = 0,
to = 1
)
approx(
d$x,
d$y,
xout = 0.95
)$y
This is dramatically faster when many calculations are needed.
36. Create a reusable Shapiro distribution object
We can take this idea one step further.
make_shapiro_distribution <- function(
n,
B = 100000,
seed = NULL
) {
if (!is.null(seed))
set.seed(seed)
W <- rshapiro(
N = B,
n = n
)
density_object <-
density(
W,
from = 0,
to = 1
)
result <- list(
n = n,
B = B,
W = W,
density = density_object
)
class(result) <-
"shapiro_distribution"
return(result)
}
Create one:
SW30 <-
make_shapiro_distribution(
n = 30,
B = 100000,
seed = 123
)
37. Use the stored distribution
Cumulative probability
mean(
SW30$W <= 0.95
)
Critical value
quantile(
SW30$W,
0.05
)
Density
approx(
SW30$density$x,
SW30$density$y,
xout = 0.95
)$y
38. Why this exercise is statistically useful
Building these functions reveals something deeper about hypothesis testing.
Every classical test requires:
- a test statistic;
- a null distribution for that statistic;
- an observed statistic;
- a tail probability.
For Shapiro-Wilk:
↓
Wobserved
↓
Compare with null distribution of W
↓
P(W ≤ Wobserved)
↓
p-value
39. The same idea appears in many statistical tests
For a t-test:
For a chi-square test:
For ANOVA:
For Shapiro-Wilk:
The only major difference is that the W distribution is more complicated and depends explicitly on sample size.
40. Final perspective
Base R does not provide:
dshapiro()
pshapiro()
qshapiro()
rshapiro()
because the Shapiro-Wilk statistic does not belong to a simple named probability family with a convenient closed-form density.
But simulation allows us to construct all four operations.
simulate W under normality
pshapiro()
estimate P(W ≤ w)
qshapiro()
find the W corresponding to a chosen probability
dshapiro()
estimate the local density of W
This exercise also exposes the machinery hidden underneath a hypothesis test.
The p-value returned by:
shapiro.test(x)
is not a mysterious number produced by software. Conceptually, it answers a very concrete question:
If the observations really came from a normal population, how often would samples of this size produce a Shapiro-Wilk statistic at least this far into the lower tail?
Monte Carlo simulation lets us watch that probability emerge directly from repeated random experiments.
And that is one of the most useful bridges between probability theory and statistical inference:
