Numbers 1-100, calculating how many I must pick to get x% of reliability I would like to know the formula or method of calculating this.
For example, I want to know how many numbers I must pick in order to have 95% confidence that I will be within 2% of the average.
I'm also interested in doing a similiar analysis for the median as well.
Update: Numbers are not put back in the pool. I think this has to do with combinations rather than trials? For example, if I pick 1 number or 99 numbers there are 100 combinations, if I pick 100 numbers there is 1 combination etc. I think you could figure out all of the combinations when you pick 1 ball, 2 balls, 3 balls etc and see how many are within 2% of the median.
Im interested in seeing just how reliable the average price statistic is for the housing market. For example, in a single month 5% of the houses may sell, or maybe 20% will. Housing prices are all over the map but I wanted to do something simplistic like this.
And now I have thought of something else...let's say I have 1-500 instead of 1-100. If I pick from 1-100 and pick 10 balls(10% of them), Im guessing that the answer would differ from taking 50 balls out of 1-500(10% of them as well)?
 A: I guess that you are looking for experimental design for estimating the mean. I am not aware of something similar for the median though.
A: I did a little simulation in R, which resulted in 98. If you go to 5% of the mean then you're down to 93 and if you go to 10% then it's 76.
Similar results are obtained for the median.
Please, read whuber's comments to my answer
The code I used:
numbers <- 1:100
numTrials <- 1000

means <- matrix(nrow = numTrials, ncol = length(numbers))

for (n in 1:length(numbers))
    {
    for (t in 1:numTrials)
        {
        vals <- sample(numbers, n)
        means[t, n] <- mean(vals) # or median if you want the median 
        }
    }

realAvg <- mean(numbers)
uplim <- realAvg + 0.05*realAvg
downlim <- realAvg - 0.05*realAvg

res <- unlist(apply(means, 2, function(x){
       length(which(x>downlim & x<uplim))/length(x)}))
plot(1:length(numbers), res*100, t="l", xlab="Numbers picked", 
       ylab="% trials with mean within limits")
abline(h=95, col="red")

print(which(res>=0.95))


PS: anyone would like to explain the "sawtooth" look of the graph?
A: A couple of points first:
1) The answer will be very large, so we are safe to assume the Central Limit Theorem is in operation, even if the distribution is all 1s and 100s (i.e. Bernoulli).
2) We can use normal deviates (z-scores) rather than t-tests because the required N will be very large - this is a simplification.
3. If we are sampling without replacement the problem statement seems to have only one interpretation - the numbers 1,...,100 - a uniform discrete distribution. whuber and nico have already provided answers, so I will deal with sampling with replacement and an undefined probability mass function (PMF).  
Approach 1:
Assume a discrete uniform distribution and sampling with replacement. The mean is 50.5, the SD is sqrt((n^2 - 1)/12) ~ 29
We want a maximum error of 2% so delta = abs(est_mean - true_mean)/true_mean = 1.01
Assume a normal distribution, the standard error of the mean is:
\begin{align}
SEM &= SD/sqrt(n)\\
sqrt(n) &= SD/SEM\\
n &= (SD/SEM)^2\\
n &= (SD/(delta/1.96))^2\\
n &= (29/(1.01/1.96))^2\\
n &\approx 3167\\
\end{align}
Approach 2:
Assume the worst case.
I'm pretty sure this is a 50:50 split in the pmf  between 1 and 100. If so our mean is 50.5 and our SD is sqrt(49.5^2)=49.5.
\begin{align}
SEM &= SD/sqrt(n)\\
sqrt(n) &= SD/SEM\\
n &= (SD/SEM)^2\\
n &= (SD/(delta/1.96))^2\\
n &= (49.5/(1.01/1.96))^2\\
n &\approx 9227\\
\end{align}
I'll let someone else do the sample sizes for medians.
A: numbers <- 1:500
numTrials <- 1000

means <- matrix(nrow = numTrials, ncol = length(numbers))

for (n in 1:length(numbers))
    {
    for (t in 1:numTrials)
        {
        vals <- sample(numbers, n)
        means[t, n] <- median(vals)
        }
    }

realAvg <- median(numbers)
uplim <- realAvg + 0.05*realAvg
downlim <- realAvg - 0.05*realAvg

res <- unlist(apply(means, 2, function(x){
       length(which(x>downlim & x<uplim))/length(x)}))
plot(1:length(numbers), res*500, t="l", xlab="Numbers picked", 
       ylab="% trials with mean within limits")
abline(h=95, col="red")

print(which(res>=0.95))

