fit.R-20170905

yinterval=as.matrix(yinterval)
f=as.matrix(yinterval[,3:18])
xint=seq(5000,155000,10000)
par(mfrow=c(2,2))
plot(xint,f[4,])

m=matrix(0,nrow(f),1)
v=matrix(0,nrow(f),1)
xmu=matrix(0,nrow(f),1)
xsigma=matrix(0,nrow(f),1)
fpdf=matrix(0,nrow(f),ncol(f))
r=matrix(0,nrow(f),ncol(f))

for (i in 1:nrow(f)){
m[i,1]=xint%*%as.numeric(f[i,])
v[i,1]=((xint-m[i,1])^2%*%as.numeric(f[i,]))^(1/2)

# xmu[i,1]=log(m[i,1]^2/(1+v[i,1]/m[i,1]^2)^(1/2))
# xsigma[i,1]=(log(1+v[i,1]/(m[i,1]^2)))^(1/2)

fpdf[i,]=dlnorm(xint,meanlog=m[i,1],sdlog=v[i,1],log=FALSE)
r[i,]=abs(fpdf[i,]-as.numeric(f[i,]))/fpdf[i,]
}

fmean=matrix(apply(f,1,mean))
rsqu[i,1]=(sum(fpdf[i,]-fmean[i]))^2/(sum(f[i,]-fmean[i]))^2

plot(xint,dlnorm(xint,m[1,1],v[1,1]),type="l",xlab="x",ylab="f(x)")
lines(density(x1),col="red")
legend("topright",c("True Density","Estimate"),lty=1,col=1:2)

par(mfrow=c(1,1))
qqplot(f[1,],fpdf[1,])
ks.test(as.numeric(f[1,]),fpdf[1,])

# x <- 1:10
# sample(c(x[x >8],x[x <8]),10,replace=TRUE,c(0.1,0.9)) 
# sample(c(1:2,3,5),100,replace = TRUE,prob =c(0.1,0.5,0.5))
# sample(x[x >8],100,replace=TRUE,prob=0.2)

xs=matrix(0,nrow(f),1e7)
for (i in nrow(f)){
 xs[i]=sample(0:10000,1e7*as.numeric(f[i,1]),replace = TRUE)
  for (j in 1:15){
    xs[i]=c(xs,sample(j*10000:(j+1)*10000,1e7*as.numeric(f[i,j+1]),replace = TRUE))
}
}

set.seed(1234)
rlnorm(1000,meanlog=m[1,1],sdlog = v[1,1])