I encourage you to follow up on my points here with a statistician in your department. Always better to get more brains looking at the same thing. There are various ways to calculate confidence intervals for series of data.
You could start with a typical formula for the standard error of the mean (SE mean), which is the standard deviation divided by the square root of the sample n:
σ / sqrt (n)
Lambda is a median value, so, we have to convert the SE mean to SE median:
1.253 * (σ / sqrt (n))
(various articles on the WWW validate this)
As you are dealing with a Chi-squared quantile distribution, in order to calculate the confidence intervals, I believe you can do this with the following:
pvalue <- runif(1000000, min=0, max=1)
chisq <- qchisq(1 - pvalue, 1)
lambda <- median(chisq) / qchisq(0.5, 1)
SE.median <- qchisq(0.975, 1) * (1.253 * ( sd(chisq) / sqrt( length(chisq) ) ) )
lambda + (SE.median / qchisq(0.5, 1))
lambda - (SE.median / qchisq(0.5, 1))
I increased the number of random p-values to 1 million just to make the confidence intervals tighter. At 1000 (your original post), they were pretty wide.
Take this as a first step into this. I encourage you to look up all functions / formulae that I have used.