\documentclass[11pt]{article}

\usepackage{indentfirst}
\usepackage{natbib}
\usepackage{url}

\begin{document}

\title{Fuzzy Error Bars}
\author{Charles J. Geyer}
\maketitle

\section{Introduction}

This document tries out ``fuzzy error bars'' using the R package \verb@ump@.
The idea was mentioned in \citet{GeyerMeeden05}, but was not adequately
explored due to lack of space.

\section{Example 1: Keeping It Simple}

For our first example, we make up some data.
<<ex1-data>>=
set.seed(42)
fred <- c("red", "orange", "yellow", "green", "blue", "indigo", "violet",
    "cyan", "magenta", "white", "black")
n <- seq(15, by = 5, length = length(fred))
x <- rbinom(length(n), n, 0.5)
@
and calculate the fuzzy confidence intervals
<<ex1-fuzzy>>=
library(ump)
p <- seq(0, 1, 0.0001)
alpha <- 0.05
phi <- matrix(NaN, length(p), length(fred))
for (i in 1:ncol(phi))
    phi[ , i] <- umpu.binom(x[i], n[i], p, alpha)
dimnames(phi) <- list(as.character(p), fred)
@

While we are at it, we calculate the core
<<ex1-core>>=
make.core <- function(phi) {
    nam <- dimnames(phi)[[2]]
    p <- as.numeric(dimnames(phi)[[1]])
    stopifnot(all(0 <= p & p <= 1))
    foo <- matrix(NaN, ncol(phi), 2)
    for (i in 1:ncol(phi)) {
        phii <- phi[, i]
        foo[i, 1] <- min(p[1 - phii == 1])
        foo[i, 2] <- max(p[1 - phii == 1])
    }
    dimnames(foo) <- list(nam, c("low", "high"))
    return(foo)
}
make.core(phi)
@
and support
<<ex1-core>>=
make.support <- function(phi) {
    nam <- dimnames(phi)[[2]]
    p <- as.numeric(dimnames(phi)[[1]])
    stopifnot(all(0 <= p & p <= 1))
    foo <- matrix(NaN, ncol(phi), 2)
    for (i in 1:ncol(phi)) {
        phii <- phi[, i]
        foo[i, 1] <- min(p[1 - phii > 0])
        foo[i, 2] <- max(p[1 - phii > 0])
    }
    dimnames(foo) <- list(nam, c("low", "high"))
    return(foo)
}
make.support(phi)
@

Figure~\ref{fig:one} is produced by the following code
<<label=fig1plot,include=FALSE>>=
par(mar = c(5, 6, 1, 1) + 0.1)
delta <- 0.3
foo <- make.support(phi)
plot(range(foo), c(0, length(fred) - delta),
    type = "n", axes = FALSE, ylab = "", xlab = "")
axis(side = 1)
box()
for (i in 1:ncol(phi)) {
    phii <- phi[ , i]
    ppp <- p[1 - phii > 0]
    xxx <- c(ppp, rev(range(ppp)))
    yyy <- (1 - phii)[1 - phii > 0]
    yyy <- c(yyy, 0, 0)
    yyy <- i - 1 + yyy * (1 - delta)
    polygon(xxx, yyy)
    lines(rep(x[i] / n[i], 2), c(i - 1, max(yyy)))
}
axis(side = 2, at = seq(along = fred) - 1 + (1 - delta) / 2, lab = fred,
    las = 1, padj = 1)
title(xlab = expression(p))
@
and appears on p.~\pageref{fig:one}.
\begin{figure}
\begin{center}
<<label=fig1,fig=TRUE,echo=FALSE>>=
<<fig1plot>>
@
\end{center}
\caption{Multiple Fuzzy Error Bars.  Fuzzy confidence intervals associated
with the UMPU test for the binomial distribution.  Sample sizes range
between \Sexpr{min(n)} and \Sexpr{max(n)}.  Height of the error bar denotes
the amount of fuzziness (full height, no fuzziness).  Line in the middle
of each bar denotes the MLE.}
\label{fig:one}
\end{figure}

\section{Example 2: Pointy Haired Bars (PHB)}

The associate editor handling \citet{GeyerMeeden05} quite rightly pointed
out that it is not \emph{always} the case that our fuzzy error bars have
full height, so the implied vertical axis for a fuzzy error bar goes from
zero at the bottom to one at the top.

However these special cases are self-labeling by their special form.
We call them \emph{pointy haired bars} (PHB), a Dilbert allusion.
To see how this works, let's do an example of the kind of confidence interval
explanation seen in many introductory statistics textbooks where a single
plot shows the result of a small simulation study, and (sure enough)
about 95\% of the 95\% confidence intervals for simulated data
cover the ``simulation truth'' parameter value and
about 5\% don't cover.

Intro stats books usually do this with normal data
(or at least quasicontinuous data).  We'll use binomial.
Also (to be fair to the intro stats books) we obey the ``more than five
expected in each cell'' rule for contingency tables with $n = 100$ and
$p_0 = 0.05$ for the simulation truth.  Thus we are in ``large $n$''
territory according to the time honored (but dumb) rule.
Some newer textbooks, such as \citet{WildSeber00} do strongly criticize
the time honored rule of thumb.

Let's do it.
<<ex2-sim>>=
n <- 100
p0 <- 0.05
nsim <- 50
x <- rbinom(nsim, n, p0)
print(x)
@
O.~K.  We got a zero, so we'll see the problem which is the ostensible
point of this section.
Figure~\ref{fig:two} is produced by the code very similar
to Figure~\ref{fig:one}, so we will not show the actual code.
<<ex2-fuzzy,echo=FALSE,hide=TRUE>>=
p <- seq(0, 1, 0.0001)
alpha <- 0.10
phi <- matrix(NaN, length(p), length(x))
for (i in 1:ncol(phi))
    phi[ , i] <- umpu.binom(x[i], n, p, alpha)
dimnames(phi) <- list(as.character(p), NULL)
@
<<label=fig2plot,include=FALSE,echo=FALSE,hide=TRUE>>=
par(mar = c(5, 1, 1, 1) + 0.1)
delta <- 0.3
foo <- make.support(phi)
plot(range(foo), c(0, length(x) - delta),
    type = "n", axes = FALSE, ylab = "", xlab = "")
axis(side = 1)
box()
for (i in 1:ncol(phi)) {
    phii <- phi[ , i]
    ppp <- p[1 - phii > 0]
    xxx <- c(ppp, rev(range(ppp)))
    yyy <- (1 - phii)[1 - phii > 0]
    yyy <- c(yyy, 0, 0)
    yyy <- i - 1 + yyy * (1 - delta)
    polygon(xxx, yyy)
    lines(rep(x[i] / n, 2), c(i - 1, max(yyy)))
}
title(xlab = expression(p))
abline(v = p0, lty = 2)
@
The figure appears on p.~\pageref{fig:two}.
\begin{figure}
\begin{center}
<<label=fig2,fig=TRUE,echo=FALSE>>=
<<fig2plot>>
@
\end{center}
\caption{Multiple Fuzzy Error Bars.  Fuzzy confidence intervals associated
with the UMPU test for the binomial distribution.
Height of the error bar denotes the amount of fuzziness
(full height, no fuzziness).  All bars with flat top have full height.
The bar without flat top
(number \Sexpr{seq(along = x)[x == 0]} counting from bottom)
only goes up to the confidence level,
which is \Sexpr{1 - alpha}.  Line in the middle
of each bar denotes the MLE.  Vertical dashed line is the ``simulation
truth'' parameter value \Sexpr{p0}.  All sample sizes $n = \Sexpr{n}$.}
\label{fig:two}
\end{figure}

Hmmmmm.  A bit cluttered with \Sexpr{nsim} intervals, but we'll leave it.
The fuzziness can still be seen.

Oh.  Almost forgot.  How well did we do?
<<ex2-eval>>=
coverage <- 1 - phi[p == 0.05, ]
sum(coverage == 1)
sum(coverage == 0)
sort(coverage[0 < coverage & coverage < 1])
mean(coverage)
sd(coverage) / sqrt(nsim)
@
So we fully covered \Sexpr{sum(coverage == 1)} times,
completely missed \Sexpr{sum(coverage == 0)} times,
and partially covered \Sexpr{sum(0 < coverage & coverage < 1)} times.
The empirical expected coverage was \Sexpr{round(mean(coverage), 4)} as opposed
to the nominal \Sexpr{1 - alpha} level of the intervals.
Not bad for simulation sample size \Sexpr{nsim}, only a little more
than one standard error (Monte Carlo standard error
\Sexpr{round(sd(coverage) / sqrt(nsim), 4)}) away from the
theoretically correct result.

\section{Example 3: Technical Quibbles}

It is tempting to describe the exceptional interval in Figure~\ref{fig:two}
as abutting the end of the data space (zero), but that's not a good
description.  Because the interval two above it, \emph{also} satisfies that
description, but is not problematic.

Too see this clearly, here's Figure~\ref{fig:two} redone to just show
the intervals for $x = 0$ or $x = 1$.
<<label=fig3plot,include=FALSE,echo=FALSE,hide=TRUE>>=
par(mar = c(5, 1, 1, 1) + 0.1)
delta <- 0.3
phi <- phi[ , x == 0 | x == 1]
x <- x[x == 0 | x == 1]
foo <- make.support(phi)
plot(range(foo), c(0, length(x) - delta),
    type = "n", axes = FALSE, ylab = "", xlab = "")
axis(side = 1)
box()
for (i in 1:ncol(phi)) {
    phii <- phi[ , i]
    ppp <- p[1 - phii > 0]
    xxx <- c(ppp, rev(range(ppp)))
    yyy <- (1 - phii)[1 - phii > 0]
    yyy <- c(yyy, 0, 0)
    yyy <- i - 1 + yyy * (1 - delta)
    polygon(xxx, yyy)
    lines(rep(x[i] / n, 2), c(i - 1, max(yyy)))
}
title(xlab = expression(p))
abline(v = p0, lty = 2)
@
This is Figure~\ref{fig:three}, which appears on p.~\pageref{fig:three}.
\begin{figure}
\begin{center}
<<label=fig3,fig=TRUE,echo=FALSE>>=
<<fig3plot>>
@
\end{center}
\caption{Fuzzy Error Bars that Abut the End.  Fuzzy confidence intervals
from Figure~\ref{fig:two} that have zero for their left endpoint.
The bar with flat top has full height.
The bar with pointed top only goes up to the confidence level,
which is \Sexpr{1 - alpha}.}
\label{fig:three}
\end{figure}

Thus we see that the correct description of \emph{full height} fuzzy
error bars is \emph{having flat top}, not \emph{abutting the end of
the range}.  There's another technical reason for this quibble as
the following case shows.
<<ex4-calc>>=
n <- 5
x <- 2
alpha <- 2 / 3
phii <- umpu.binom(x, n, p, alpha)
@
<<label=fig4plot,include=FALSE,echo=FALSE,hide=TRUE>>=
par(mar = c(5, 1, 1, 1) + 0.1)
delta <- 0.3
phi <- as.matrix(phii)
dimnames(phi) <- list(as.character(p), NULL)
foo <- make.support(phi)
plot(range(foo), c(0, length(x) - delta),
    type = "n", axes = FALSE, ylab = "", xlab = "")
axis(side = 1)
box()
for (i in 1:ncol(phi)) {
    phii <- phi[ , i]
    ppp <- p[1 - phii > 0]
    xxx <- c(ppp, rev(range(ppp)))
    yyy <- (1 - phii)[1 - phii > 0]
    yyy <- c(yyy, 0, 0)
    yyy <- i - 1 + yyy * (1 - delta)
    polygon(xxx, yyy)
    lines(rep(x[i] / n, 2), c(i - 1, max(yyy)))
}
title(xlab = expression(p))
@
This is Figure~\ref{fig:four}, which appears on p.~\pageref{fig:four}.
\begin{figure}
\begin{center}
<<label=fig4,fig=TRUE,echo=FALSE>>=
<<fig4plot>>
@
\end{center}
\caption{Another Fuzzy Error Bar without Full Height.
Fuzzy confidence interval associated
with the UMPU test for the binomial distribution.
Data $x = \Sexpr{x}$, $n = \Sexpr{n}$, confidence level
\Sexpr{round(100 * (1 - alpha), 2)}\%.
The value of the membership function of the fuzzy interval
goes from zero at the bottom to \Sexpr{round(max(1 - phii), 3)} at the
apex.}
\label{fig:four}
\end{figure}

Admittedly a
\Sexpr{round(100 * (1 - alpha), 2)}\% confidence interval
has no conceivable application.  We show it only to make the technical
point that very low confidence level combined with small $n$ can give
another case where the fuzzy error bar does not have full height.
And it is also another case where \emph{full height} goes with \emph{flat top}
and \emph{non-full height} or \emph{empty core} goes with \emph{pointed top}.

For the curious, the reason for the pointed top is that in order to reject
at such a high $\alpha$ we must reject, at least partially, at \emph{all}
points, even the point $x$, even when $p = \hat{p}$, where the fuzzy
confidence interval membership function always achieves its maximum
(possibly one, possibly less than one).  This latter claim is a theorem from
unpublished notes of the author: the membership function of a fuzzy confidence
interval associated with the UMPU test for a one-parameter discrete exponential
family achieves its maximum when the parameter value equals the MLE.

After that bit of theory, we see that abstractly the behavior in
the bottom interval of Figure~\ref{fig:three} is just the same as
the behavior of the interval in Figure~\ref{fig:four}.  The only
difference is that in Figure~\ref{fig:three}, the MLE of the mean value
parameter is on the boundary
of the sample space, and when that happens, the MLE distribution
concentrates all probability at the point on the boundary, and thus
we see the irregularity for all $\alpha$ levels.  But when the
phenomenon occurs in the interior of the sample sample space, we see
the irregularity only for high $\alpha$ levels, hence low confidence levels
$$
   1 - \alpha < f_{\hat{\mu}}(x)
$$
where $\hat{\mu} = x$ is the MLE of the mean value parameter.
See equations (3.7a) and (3.7b) and the following discussion
in \citet{GeyerMeeden05} for explanation.

\begin{thebibliography}{}

\bibitem[{Geyer and Meeden(2005)}]{GeyerMeeden05}
\textsc{Geyer, C.~J.} and \textsc{Meeden, G.~D.} (2005).
\newblock Fuzzy and randomized confidence intervals and p-values.
\newblock Revised and resubmitted (thrice) to \emph{Statistical Science},
  \url{http://www.stat.umn.edu/geyer/fuzz}.

\bibitem[{Wild and Seber(2000)}]{WildSeber00}
\textsc{Wild, C.~J.} and \textsc{Seber, G.~A.~F.} (2000).
\newblock \emph{Chance Encounters: A First Course in Data Analysis
    and Inference}.
\newblock Wiley.

\end{thebibliography}

\end{document}