We simulate the sampling distribution of the average and the median in a sample generated from the Normal distribution. In order to do so we copy the code that was used in Subsection Subsection 10.3.2, replacing the object “mid.range" by the object “X.med" and using the function “median" in order to compute the sample median instead of the computation of the mid-range statistic:
mu <- 3
sig <- sqrt(2)
X.bar <- rep(0,10^5)
X.med <- rep(0,10^5)
for(i in 1:10^5) {
X <- rnorm(100,mu,sig)
X.bar[i] <- mean(X)
X.med[i] <- median(X)
}
The sequence “X.bar" represents the sampling distribution of the sample average and the sequence “X.med" represents the sampling distribution of the sample median. We apply the function “mean" to these sequences in order to obtain the expectations of the estimators:
The expectation of the measurement, the parameter of interest is equal to 3. Observe that expectations of the estimators are essentially equal to the expectation[^13]. Consequently, both estimators are unbiased estimators of the expectation of the measurement.
Observe that the variance of the sample average is essentially equal to \(0.020\) and the variance of the sample median is essentially equal to \(0.0312\text{.}\) The mean square error of an unbiased estimator is equal to its variance. Hence, these numbers represent the mean square errors of the estimators. It follows that the mean square error of the sample average is less than the mean square error of the sample median in the estimation of the expectation of a Normal measurement.
We repeat the same steps as before for the Uniform distribution. Notice that we use the parameters \(a=0.5\) and \(b=5.5\) the same way we did in Subsection Subsection 10.3.2. These parameters produce an expectation \(\Expec(X) = 3\) and a variance \(\Var(X) = 2.083333\text{:}\)
a <- 0.5
b <- 5.5
X.bar <- rep(0,10^5)
X.med <- rep(0,10^5)
for(i in 1:10^5) {
X <- runif(100,a,b)
X.bar[i] <- mean(X)
X.med[i] <- median(X)
}
Applying the function “mean" to the sequences that represent the sampling distribution of the estimators we obtain that both estimators are essentially unbiased[^14]:
Observe \(0.021\) is, essentially, the value of the variance of the sample average[^15]. The variance of the sample median is essentially equal to 0.061. The variance of each of the estimators is equal to it’s mean square error. This is the case since the estimators are unbiased. Consequently, we again obtain that the mean square error of the sample average is less than that of the sample median.
Assuming that the file “ex2.csv" is saved in the working directory, one may read the content of the file into a data frame and produce a summary of the content of the data frame using the code:
## id sex age bmi
## Min. :1024982 FEMALE:74 Min. :26.00 Min. :15.12
## 1st Qu.:3172782 MALE :76 1st Qu.:32.00 1st Qu.:22.02
## Median :5200484 Median :35.00 Median :25.16
## Mean :5463304 Mean :35.09 Mean :24.76
## 3rd Qu.:7982902 3rd Qu.:37.00 3rd Qu.:27.49
## Max. :9934175 Max. :45.00 Max. :35.24
## systolic diastolic group
## Min. :100.8 Min. : 51.98 HIGH : 37
## 1st Qu.:118.1 1st Qu.: 75.02 LOW : 3
## Median :124.3 Median : 83.19 NORMAL:110
## Mean :125.3 Mean : 82.44
## 3rd Qu.:132.6 3rd Qu.: 88.83
## Max. :154.5 Max. :112.71
Examine the variable “group". Observe that the sample contains 37 subjects with high levels of blood pressure. Dividing 37 by the sample size we get:
Notice that the expression “ex2$group == HIGH" produces a sequence of length 150 with logical entries. The entry is equal to “TRUE" if the equality holds and it is equal to “FALSE" if it dos not[^17]. When the function “mean" is applied to a sequence with logical entries it produces the relative frequeny of the TRUEs in the sequence. This corresponds, in the corrent context, to the sample proportion of the level “HIGH" in the variable “ex2$group".
Make sure that the file “pop2.csv" is saved in the working directory. In order to compute the proportion in the population we read the content of the file into a data frame and compute the relative frequency of the level “HIGH" in the variable “group":
The simulation of the sampling distribution involves a selection of a random sample of size 150 from the population and the computation of the proportion of the level “HIGH" in that sample. This procedure is iterated 100,000 times in order to produce an approximation of the distribution:
P.hat <- rep(0,10^5)
for(i in 1:10^5) {
X <- sample(pop2$group,150)
P.hat[i] <- mean(X == "HIGH")
}
mean(P.hat)
## [1] 0.2811699
Observe that the sampling distribution is stored in the object “P.hat". The function “sample" is used in order to sample 150 observation from the sequence “pop2$group". The sample is stored in the object “X". The expression “mean(X == HIGH)" computes the relative frequency of the level “HIGH" in the sequence “X".
At the last line, after the production of the sequence “P.hat" is completed, the function “mean" is applied to the sequence. The result is the expected value of estimator \(\hat P\text{,}\) which is equal to \(0.2812307\text{.}\) This expectation is essentially equal to the probability of the event \(p = 0.28126\text{.}\)[^18]
We get that the proposed variance in Section 5 is \(0.0013476850\text{,}\) which is in good agreement with the value \(0.00135\) that was obtained in the simulation[^19].
