Horseshoes, Hand grenades and Slow-dancing

Things where close counts!

Hi All. This morning I had an idea about a distribution for data where close was good enough. For example what if you are trying to model where users touch on a screen given that any place they touch a target is equivalent but any where off of the target is bad. There is no reason for the the distribution to be peaked like a Gaussian; so estimates from peaked distributions will have incorrect variances. Also It is nice for the function to be made of simple functions so that it is not computationally intensive to calculate the cdf / pdf.

The basic idea for distribution came from a Prof I had in my undergrad who used the presented cdf (with ) to plot limits for students. Note I was hoping to find a better web page to link to  but Google is obsessed with some yob that is a professor at Harvard with the same name. So I call this the “Close-enough” distribution. Enough rambling!

The CDF is a piece wise function where if in (, ),  if , and if .

The cumulative density function for Close-enough distribution.

The pdf is where in (, ),  other wise.

The pdf of the Close-enough distribution.

I would like to like to provide a most likely estimate of however I am lazy and MCMC is a good way to get an estimate.Code for the MCMC done here is given below.

The PPD of estimates for a simulated example of data drawn from the Close-enough distribution with true scale parameter given by the vertical line.

Code:

g.pdf <-function (a, x=seq(-11,11, length.out=10001) )
{
f1 = 1/(4*a) + 0*x
f2 = (a/4)*x^(-2)
f3 = (a/4)*x^(-2)

f = f1
f[x>a] = f2[x>a]
f[x<(-a)] = f3[x<(-a)]

#plot(c(-11,x,11), c(0,f,0), type=”n”, xlim=c(-10, 10), xaxs=”i”, xlab=”x”, ylab=”P(x)”  )
#polygon(c(-11,x,11), c(0,f,0), col=”grey”)

return(f)

}

g.inv.cdf <-function (a, P)
{
tv = 4*a*P – 2*a

if ((tv >= -a ) & (tv <= a ))
{
x = tv
}

if ((tv > a ))
{
x = -a^2/(4*a*(P-1))
}

if ((tv < -a ))
{
x = -a/(4*P)
}

return(x)
g.rced <-function (n= 1000, A=1)
{
out = 1:n

for (i in 1:n)
{
out[i] = g.inv.cdf(P= runif(1), a=A)

}
#out = out[out>-10]
#out = out[out<10]

h=hist(out, breaks=101, freq=F)
d= g.pdf(x=h$mids, a=A)

lines(h$mids, d, col=”blue”)
return(out)
}

g.mcmc <-function (N = 100, n.=100)
{
d = g.rced(n=n., A=1)

LL= -Inf
a = 1

mod  = c(LL, a)
rec  = data.frame(LL= 1:N, a= 1:N)
Q.sd = 0.1

for (i in 1:N )
{
print(i)
mod.new= mod
mod.new[2] = mod.new[2] + rnorm(1, sd= Q.sd, mean=0)

mod.new[1] = sum(log(g.pdf(x=d, a = mod.new[2])))
A = exp(mod.new[1]-mod[1])
r = runif(1)

if (r<=A)
{
mod = mod.new
}
rec[i,c(1,2)] = mod
}
#print(rec)
h=hist(rec[,2], breaks=101, plot = F)

plot(h$mids, h$dens, xlab=”a”, ylab=”P(a)”,type=”n”)
polygon(c(h$mids[1], h$mids, h$mids[2]), c(0,h$dens,0), col=”grey”, border=F)
abline(v=1, lwd=2, col=”blue”)

}