#(b = (mean(x)^2)/(var(x)+mean(x)^2))
#k = seq(0.00001,10,length=10000)
#gk = exp(lgamma(1+1/k)*2) - b*exp(lgamma(1+2/k))
#plot(gk~k,type='l')
#abline(h=0, col=2,lwd=.1)
#sele = !is.na(gk)
#gk2 = abs(gk[sele])
#k2 = k[sele]
#alpha = k[which(gk2==min(gk2))]
#lambda = mean(x)/gamma(1+1/alpha)
#c(alpha,lambda)
data(birth)
(ofc = binning(birth$Head))
(bwt = binning(birth$Weight, scale=100))
(out0 = bfmm(ofc, m=1))
mu = 34.5; s = 1.5
y = rnorm(100,mu,s) #raw data
x = round(y) #rounded data
binx = binning(x)
x0 = seq(min(y)-sd(y),max(y)+sd(y),length=200)
f0 = dnorm(x0,mu,s)
## Fit a single well-known distribution to binned data
xmin=29; xmax = 41;
hist(binx,xlim=c(29,41))
lines(f0~x0, col = 1)
lines(density(x+runif(length(x))-.5), col=1,lty=2)
out1 = bfmm(x, m=1, type='normal')
lines(density(out1), col=2)
out2 = bfmm(x, m=1, type='gamma')
lines(density(out2), col=4)
out3 = bfmm(x, m=1, type='weibull')
lines(density(out3), col=3)
out4 = bfmm(x, m=1, type='beta')
lines(density(out4), col=5)
out5 = bfmm(x,m=2)
lines(density(out5), col=2,lwd=3)
legend(xmax, max(binx$counts)/sum(binx$counts),
xjust=1,yjust=1,
legend=c("true","kde(x+u)","normal","gamma","weibull","beta","normal m=2"),
col = c(1,1,2,4,3,5,2), lty=c(1,2,1,1,1,1,1),lwd=c(1,1,1,1,1,1,3))
## Fit with different methods
hist(binx,xlim=c(29,41))
lines(f0~x0, col = 1)
lines(density(x+runif(length(x))-.5), col=1,lty=2)
out6 = bfmm(x, m=2, method='nelder')
lines(density(out6), col=2,lwd=1)
out7 = bfmm(x, m=2, method='newton')
lines(density(out7), col=3,lwd=1)
out8 = bfmm(x, m=1, method='nelder')
lines(density(out8), col=4,lwd=1)
out9 = bfmm(x, m=1, method='newton')
lines(density(out9), col=5,lwd=1)
out10 = bfmm(x, mu=c(34,36), method='newton')
lines(density(out10), col=6,lwd=1)
Run the code above in your browser using DataLab