Post date: Jul 07, 2015 7:42:10 AM
Is it possible to derive estimates by using the entire empirical distribution function (EDF) as the statistic rather than focusing on moments or other indexing characterizations? The calculations below with the normal distribution illustrate that it is possible to do this. However, it is not clear that there is any advantage in using this approach. The calculations include an estimate based on the whole EDF and an estimate based on the sample mean and sample standard deviation. The former seems to be uniformly tighter than the latter. Missing from these example calculations is the analogous two-dimensional plausibility map corresponding to the c-boxes based on the chi squared and t-distributions. <<Perhaps Michael can advise what these should look like, if they exist.>>
The code shown below takes little time to run and it crudely shows the differences between estimating with moments and estimating with the whole EDF. To see finer detail, use instead the code attached to the bottom of this page, which memoises some results previously computed.
# fit a normal with the entire empirical distribution function, i.e., the original sample data
cbox.normal <- function(
x, # data
mus = -20:20 * sd(x) / 10 + mean(x), # possible parameter values we're considering
sigmas = seq(1/4, 4,length=20) * sd(x), # possible parameter values we're considering
many = 1000, # number of bootstrap replications
steps = 1000 # steps in the distributional approximations
) {
areametric = function(cdf, x) {
# maybe we should be comparing quasi-inverses rather than the cdfs
x = sort(x)
edf = rep(0,length(z))
for (k in 1:s) edf = ifelse(z<=x[[k]],edf,edf+1)
edf = edf/max(edf)
#plot(z,edf,type='S')
#rug(x)
#lines(rep(max(x),2),c(0,1))
sum(abs(cdf-edf)) # maybe should also add the tail areas below min(x) and above max(x) for this mu and sigma
}
s = length(x)
w = max(x)-min(x)
z = seq(min(x)-w,max(x)+w,length.out=steps)
ii = 1:length(mus)
jj = 1:length(sigmas)
pls = matrix(0,nrow=length(ii),ncol=length(jj))
for (j in jj) for (i in ii) {
mu = mus[[i]]
sigma = sigmas[[j]]
cdf = pnorm(z,mu,sigma)
dobserved = areametric(cdf,x) # difference between theoretical and observed distributions
dsampling = NULL
for (k in 1:many) {
y = rnorm(s,mu,sigma)
dsampling = c(dsampling, areametric(cdf, y))
}
pls[i,j] = length(dsampling[dobserved <= dsampling])
}
# normalize the plausibilities
pls = pls / many
if (max(pls)<1) cat('Plausibilities do not achieve unity\n')
list(mu=mus, sigma=sigmas, pls=pls / max(pls))
}
# now fit a normal with reference to the sample mean and sample standard deviation
cbox.Normal <- function( # notice the capital N in the function name
x, # data
mus = -20:20 * sd(x) / 10 + mean(x), # possible parameter values we're considering
sigmas = seq(1/4, 4,length=20) * sd(x), # possible parameter values we're considering
many = 1000, # number of bootstrap replications
steps = 1000 # steps in the distributional approximations
) {
metric = function(m,s, x) return(sqrt((m-mean(x))^2 + (s-sd(x))^2))
s = length(x)
w = max(x)-min(x)
z = seq(min(x)-w,max(x)+w,length.out=steps)
ii = 1:length(mus)
jj = 1:length(sigmas)
pls = matrix(0,nrow=length(ii),ncol=length(jj))
for (j in jj) for (i in ii) {
mu = mus[[i]]
sigma = sigmas[[j]]
cdf = pnorm(z,mu,sigma)
dobserved = areametric(mu,sigma,x) # difference between theoretical and observed distributions
dsampling = NULL
for (k in 1:many) {
y = rnorm(s,mu,sigma)
dsampling = c(dsampling, metric(mu,sigma, y))
}
pls[i,j] = length(dsampling[dobserved <= dsampling])
}
# normalize the plausibilities
pls = pls / many
if (max(pls)<1) cat('Plausibilities do not achieve unity\n')
list(mu=mus, sigma=sigmas, pls=pls / max(pls))
}
# this is the likelihood function traditionally used in this estimation
lkhd.normal <- function(x, mus = -20:20*sd(x)/10+mean(x), sigmas=seq(1/4,4,length=20)*sd(x)) {
ii = 1:length(mus)
jj = 1:length(sigmas)
pls = matrix(1,nrow=length(ii),ncol=length(jj))
for (j in jj) for (i in ii) {
mu = mus[[i]]
sigma = sigmas[[j]]
pls[i,j] = prod(dnorm(x, mu, sigma))
}
# normalize the plausibilities
pls = pls / many
list(mu=mus, sigma=sigmas, pls=pls / max(pls))
}
cb = cbox.normal(x) # c-box on the EDF
CB = cbox.Normal(x) # c-box on moments
safecb = cb
safeCB = CB
dput(x)
dput(cb)
dput(CB)
go = function(phi = 20,theta = 55, W=2, who=4, ...) {
if (who %in% c(1,4))
persp(cb$mu,cb$sigma,cb$pls,tickt='detailed',phi=phi,theta=theta,xlab='Normal mean',ylab='Normal stdev',zlab='Plausibility', col='lightblue',shade=0.65,ltheta=20,...)->res
Sys.sleep(W)
if (who %in% c(2,4))
persp(CB$mu,CB$sigma,CB$pls,tickt='detailed',phi=phi,theta=theta,xlab='Normal mean',ylab='Normal stdev',zlab='Plausibility', col='lightgreen',shade=0.65,ltheta=20,...)->res
Sys.sleep(W)
}
par(mfrow=c(1,1))
go()
for (phi in 0:360) go(phi,55,0)
for (t in 0:360) go(10,t-120,0,who=3)
go(5,80,4)
par(ask=TRUE)
# CLICK ON PLOT TO CONTINUE
p=5;t=80;persp(cb$mu,cb$sigma,cb$pls-CB$pls,tickt='detailed',phi=p,theta=t,xlab='Normal mean',ylab='Normal stdev',zlab='Plausibility',col='lightblue',shade=0.65,ltheta=20)->res
par(mfrow=c(3,1))
xlim=c(5.5,6.5)
ylim=c(0.5,1.5)
Z=cb$pls;contour(cb$mu,cb$sigma,Z,xlab='Normal mean',ylab='Normal stdev',xlim=xlim,ylim=ylim,col='blue')
Z=CB$pls;contour(cb$mu,cb$sigma,Z,xlab='Normal mean',ylab='Normal stdev',xlim=xlim,ylim=ylim,col='red')
Z=cb$pls;contour(cb$mu,cb$sigma,Z,xlab='Normal mean',ylab='Normal stdev',xlim=xlim,ylim=ylim,col='blue')
Z=CB$pls;contour(cb$mu,cb$sigma,Z,xlab='Normal mean',ylab='Normal stdev',xlim=xlim,ylim=ylim,col='red',add=TRUE)