Expected number of trials for all 195 numbers to be chosen in the range of 1 to 195 - Page 2 - Gambling and Probability

Two Plus Two Forums Other Topics Probability

Expected number of trials for all 195 numbers to be chosen in the range of 1 to 195

Post Reply Subscribe

...

Page 2 of 4

1 2 3 4

Page 2 of 4

1 2 3 4

12-16-2012 , 03:55 PM

#26

Sherman

Carpal \'Tunnel

Join Date: Jun 2005 Posts: 7,766

Quote:

Originally Posted by BruceZ

It should work then. What lower and upper endpoints did you use? Don't make the lower endpoint less than 195 because you can't get 195 countries in less than 195 days. If you make it too small it fails.

Yeah, I literally just copied and pasted your code. Ran it. Got an error.

Code:

> prob = function(N) {
+   s1 <- seq(1,NCoupons-1)
+   x <- sum((-1)^(s1+1) * choose(NCoupons,s1)*((NCoupons-s1)/NCoupons)^N)
+   x - 0.5
+ }
> 
> NCoupons = 195
> uniroot(prob, lower = 195, upper = 2000, tol = 0.5)
Error in uniroot(prob, lower = 195, upper = 2000, tol = 0.5) : 
  f() values at end points not of opposite sign

But if I change the lower value to something (like 300) it works:

Code:

> prob = function(N) {
+   s1 <- seq(1,NCoupons-1)
+   x <- sum((-1)^(s1+1) * choose(NCoupons,s1)*((NCoupons-s1)/NCoupons)^N)
+   x - 0.5
+ }
> 
> NCoupons = 195
> uniroot(prob, lower = 300, upper = 2000, tol = 0.5)
$root
[1] 1099.168

$f.root
[1] -2.717433e-06

$iter
[1] 6

$estim.prec
[1] 0.2544312

Quote

12-16-2012 , 04:04 PM

#27

BruceZ

Carpal \'Tunnel

Join Date: Sep 2002 Posts: 11,877

I don't know what to tell you. I restarted R, then copied and pasted exactly what you said wasn't working into the interpreter, and it worked.

> prob = function(N) {
+ s1 <- seq(1,NCoupons-1)
+ x <- sum((-1)^(s1+1) * choose(NCoupons,s1)*((NCoupons-s1)/NCoupons)^N)
+ x - 0.5
+ }
>
> NCoupons = 195
> uniroot(prob, lower = 195, upper = 2000, tol = 0.5)
$root
[1] 1099.173

$f.root
[1] -1.086849e-05

$iter
[1] 5

$estim.prec
[1] 0.25

Quote

12-17-2012 , 04:25 AM

#28

David Sklansky

Administrator

Join Date: Aug 2002 Posts: 17,352

You can estimate the answer this way:

A particular country is even money to be chosen after about 135 days. The chances that a country will be oicked after 135x8=1080 tries is 255/256. 70% of 255 is 178. So in 1080 tries its a little better than even money that one of the 195 countries wil be unpicked.

Quote

12-17-2012 , 08:00 PM

#29

BruceZ

Carpal \'Tunnel

Join Date: Sep 2002 Posts: 11,877

Disregard previous question. I got it straightened out and revised this yet again.

Quote:

Originally Posted by David Sklansky

1080 days gives a probability of 53.5% that some country will remain unpicked. The actual median is 1099. If we had 178 countries, the median would be 987, and your method gives 984 (not what I said it gave before).

David's 70% comes from ln(2) which is about 0.693. We know that the probability of not getting a particular country in N days is exactly (194/195)^N. That is about 1/2 when N is 135. It is 1/256 when N is 8*135 = 1080. That's still for 1 country, so if we treat all the countries as independent, (255/256)^N = 0.5 when N = 177. He got 178 by multiplying 256 by 0.7,and that gives 49.8%. The 177 actually comes from multiplying 255 by ln(2), or 256 by ln(2). Anyway, the basis for this estimate is

(255/256)^N = 0.5

(256/255)^N = 2

N = ln(2) / ln(256/255)

= ln(2) / ln(1 + 1/255)

=~ ln(2) / (1/255)

=~ ln(2) * 255.

=~ 0.693 * 255 =~ 176.7 or 177

=~ 0.7 * 255 = 178.5

Since we have 195 countries and not 178 countries, it will be a little higher than 50% to not have all countries in 1080 days, and it is in fact 53.5%.

Note that we could have estimated the median as 1103 (off by just 4) by solving for N in just this 3 term formula:

195*(194/195)^N - C(195,2)*(193/195)^N + C(195,3)*(192/195)^N = 0.5

But David's estimate doesn't require a calculator.

I had some stuff here before about how this is similar to the way I estimate the median for streaks without a streak calculator, and I had a demonstration of how accurate that is, and some R scripts to do that. But David's is a little different because he is applying the method a second time over the countries, plus you can estimate it with nothing but multiplications. For the streak calculation the probabilities are usually too big to approximate the log that way. What I said about the streak calculation wasn't anything I hadn't posted before, so it's gone now.

Last edited by BruceZ; 12-17-2012 at 08:49 PM.

Quote

12-18-2012 , 01:38 PM

#30

BruceZ

Carpal \'Tunnel

Join Date: Sep 2002 Posts: 11,877

We can adapt David's method to estimate the median to within 2 at 1097 without solving a polynomial. You just need a calculator that can do logs.

1. Compute 195 / ln(2) =~ 281.3

2. Compute ln(281.3) / ln(2) =~ 8.14

3. Compute ln(0.5) / ln(194/195) =~ 134.8

4. Multiply 8.14 * 134.8 =~ 1097

So in relation to David's, we found a number that about 69.3% of was 195 instead of 178, and that was 281.3 instead of 256. Then we found the power of 2 that gave that number which was 8.14 instead of 8. Then we multiplied that by 134.8 instead of 135.

Last edited by BruceZ; 12-18-2012 at 02:28 PM.

Quote

12-20-2012 , 12:17 AM

#31

BruceZ

Carpal \'Tunnel

Join Date: Sep 2002 Posts: 11,877

It is most likely that your problems with uniroot were not problems with uniroot, but with numerical instability of inclusion-exclusion. I found conditions under which inclusion-exclusion will become numerically unstable and return bogus values or overflow/underflows. One condition is when NCoupons is larger at 2838. This can be fixed by summing only 100 terms in the function instead of NCoupon terms. Even then it still becomes unstable if the lower endpoint is too small. The median for 2838 is 23604, and it worked with a lower limit down to 13000, but failed at 12000. The numerical instability is caused by subtracting large numbers that are almost equal and losing the result to finite precision. It's possible that there are islands of stability and instability rather than failing at a particular number, though I haven't confirmed this. Using interval values not ridiculously far from the right value should work ok. It should be possible to trap error values in code. The only way I could see it working for me and not for you is if there is a random component to the points that uniroot checks, or if you have a different version of uniroot.

If you want to pass a parameter like NCoupons to uniroot, you put it in the call to uniroot, not in the call to the function.

Code:

> prob = function(N, NCoupons) {
+   s1 <- seq(1,100)
+   x <- sum((-1)^(s1+1) * choose(NCoupons,s1)*((NCoupons-s1)/NCoupons)^N)
+   x - 0.5
+ }
> 
> uniroot(prob, NCoupons = 195, lower = 195, upper = 2000, tol = 0.5)
$root
[1] 1099.173

$f.root
[1] -1.086849e-05

$iter
[1] 5

$estim.prec
[1] 0.25

It passes values from the interval from lower to upper to the first parameter N, and then passes the other parameters you specify. You can also specify the interval as a vector like c(195,2000) as the 2nd argument to uniroot instead of specifying the upper and lower parameters.

Last edited by BruceZ; 12-20-2012 at 01:01 AM.

Quote

12-20-2012 , 12:37 AM

#32

Sherman

Carpal \'Tunnel

Join Date: Jun 2005 Posts: 7,766

Thanks for looking into that Bruce. I had seen something (after posting before) about passing other parameters to uniroot() so that makes sense. It is still annoying that it is unstable, but I can at least understand why now. It sounds to me like the key is to fiind a way to select intervals that are not far away...maybe using one of shortcuts you and David mentioned above would be a good place to start? Or perhaps if there was a way to identify reasonable upper and lower bounds.

To be clear, it would be nice to add this subfunction (your prob and uniroot call) inside the coupon() function that I already wrote so that you would get the Mean, Median, SD, and Variance all in one function. But to do that, you would want to safeguard yourself from this problem. One obvious choice is to have the user input values when he/she desires the median as well, but that isn't very fancy now is it? Perhaps reasonable values could be selected by using the +1 and -1 SD from the mean...I'm not sure how often (if ever) the median would be that far from the mean.

edit: Maybe something like this:

Code:

coupon2 <- function(NCoupons, MEDIAN=TRUE, tol=.001) {
  p_i <- (NCoupons - (seq(1,NCoupons) - 1)) / NCoupons
  EV <- sum(1/p_i)
  VAR <- sum((1-p_i)/p_i^2)
  SD <- sqrt(VAR)  
  if(MEDIAN==TRUE) {
    med.fun <- function(N, NCoupons) {
      s1 <- seq(1,NCoupons-1)
      x <- sum((-1)^(s1+1) * choose(NCoupons,s1)*((NCoupons-s1)/NCoupons)^N)
      x - .5
    }
    MED <- uniroot(med.fun, NCoupons=NCoupons, c(EV-SD, EV+SD), tol=tol)$root
    res <- cbind(EV, MED, VAR, SD)
    rownames(res) <- c("Results")
    colnames(res) <- c("Mean", "MED", "VAR", "SD")
    return(res)  
  }
  res <- cbind(EV, VAR, SD)
  rownames(res) <- c("Results")
  colnames(res) <- c("Mean", "VAR", "SD")
  return(res)  
}

Kinda nice.

Code:

> coupon2(195)
            Mean      MED      VAR       SD
Results 1141.292 1099.167 61212.83 247.4123

Last edited by Sherman; 12-20-2012 at 12:53 AM.

Quote

12-20-2012 , 01:34 AM

#33

BruceZ

Carpal \'Tunnel

Join Date: Sep 2002 Posts: 11,877

My version of the DS method with logs is basically dead nuts on, perhaps off by 2 or 3. The function is below. C is the number of coupons. The name dsc is for DS with a calculator. The other one ds is DS without a calculator. That should return lower and upper bounds that you could use. Reduce NCoupons-1 to something reasonable in med.fun. 100 terms should be way more than enough.

Code:

dsc = function(C) {
  log(C/log(2))/log(2) * log(0.5)/log((C-1)/C)
}

ds = function(C) {
    lower = ceiling(log(C)/log(2)) * floor(0.7*(C-1))
    upper = (ceiling(log(C)/log(2)) + 1) * floor(0.7*(C-1))
    c(lower,upper)  
}

Quote

12-20-2012 , 01:43 AM

#34

Sherman

Carpal \'Tunnel

Join Date: Jun 2005 Posts: 7,766

Cool. So if I understand correctly, you are saying that rather than using N-1 coupons, I should just change the value to a constant like 100 (when making s1)?

This is what I have now:

Code:

coupon2 <- function(NCoupons, MEDIAN=TRUE, tol=.001) {
  p_i <- (NCoupons - (seq(1,NCoupons) - 1)) / NCoupons
  EV <- sum(1/p_i)
  VAR <- sum((1-p_i)/p_i^2)
  SD <- sqrt(VAR)  
  if(MEDIAN==TRUE) {
    med.fun <- function(N, NCoupons) {
      seq.max <- ifelse(NCoupons > 100, 100, NCoupons)
      s1 <- seq(1,seq.max)
      x <- sum((-1)^(s1+1) * choose(NCoupons,s1)*((NCoupons-s1)/NCoupons)^N)
      x - .5
    }
    lower <- ceiling(log(NCoupons)/log(2)) * floor(0.7*(NCoupons-1))
    upper <- (ceiling(log(NCoupons)/log(2)) + 1) * floor(0.7*(NCoupons-1))
    MED <- uniroot(med.fun, NCoupons=NCoupons, c(lower, upper), tol=tol)$root
    res <- cbind(EV, MED, VAR, SD)
    rownames(res) <- c("Results")
    colnames(res) <- c("Mean", "MED", "VAR", "SD")
    return(res)  
  }
  res <- cbind(EV, VAR, SD)
  rownames(res) <- c("Results")
  colnames(res) <- c("Mean", "VAR", "SD")
  return(res)  
}

edit: Had to make a minor change to the code above. An error was produced when NCoupons < 100. That is fixed by the seq.max <- ifelse() part now.

Quote

12-20-2012 , 01:47 AM

#35

BruceZ

Carpal \'Tunnel

Join Date: Sep 2002 Posts: 11,877

Right, even when I did 2838 coupons, I only needed less than 10 terms.

Quote

12-20-2012 , 01:51 AM

#36

Sherman

Carpal \'Tunnel

Join Date: Jun 2005 Posts: 7,766

Quote:

Originally Posted by BruceZ

Right, even when I did 2838 coupons, I only needed less than 10 terms.

Ah...mine craps out at 2049. But seems to work (for the median) at all values below that.

Quote

12-20-2012 , 02:00 AM

#37

Sherman

Carpal \'Tunnel

Join Date: Jun 2005 Posts: 7,766

Actually, if I use the lame +1SD and -1SD for the upper and lower bounds, mine works for 2838 (and 10,000 for that matter). So I guess I am going to still with that method.

Code:

coupon3 <- function(NCoupons, MEDIAN=TRUE, tol=.001) {
  p_i <- (NCoupons - (seq(1,NCoupons) - 1)) / NCoupons
  EV <- sum(1/p_i)
  VAR <- sum((1-p_i)/p_i^2)
  SD <- sqrt(VAR)  
  if(MEDIAN==TRUE) {
    med.fun <- function(N, NCoupons) {
      seq.max <- ifelse(NCoupons > 100, 100, NCoupons)
      s1 <- seq(1,seq.max)
      x <- sum((-1)^(s1+1) * choose(NCoupons,s1)*((NCoupons-s1)/NCoupons)^N)
      x - .5
    }
    MED <- uniroot(med.fun, NCoupons=NCoupons, c(EV-SD, EV+SD), tol=tol)$root
    res <- cbind(EV, MED, VAR, SD)
    rownames(res) <- c("Results")
    colnames(res) <- c("Mean", "MED", "VAR", "SD")
    return(res)  
  }
  res <- cbind(EV, VAR, SD)
  rownames(res) <- c("Results")
  colnames(res) <- c("Mean", "VAR", "SD")
  return(res)  
}

Code:

> coupon3(10000)
            Mean      MED       VAR       SD
Results 97876.06 95767.41 164385531 12821.29

Quote

12-20-2012 , 02:06 AM

#38

BruceZ

Carpal \'Tunnel

Join Date: Sep 2002 Posts: 11,877

Wow, I just did 2838 coupons, and the lower bound from ds of 23820 FAILED, but making it even LOWER to 23000 worked. One gives a value of -2.795797e+274
and the other gives -5.451751e+272. Even with 10 terms it failed. This is a real problem if it fails that tight. It could fail anywhere. Maybe even if you use dsc for the bounds. Of course you could just return dsc in that case and it will be almost exact.

Quote

12-20-2012 , 02:28 AM

#39

Sherman

Carpal \'Tunnel

Join Date: Jun 2005 Posts: 7,766

Is the dsc number always below the exact value? It seems like it always is (for the numbers I have tried so far). Didn't try an actual test for thousands of numbers though.

Could be used as the lower bound if that was the case. Not sure what I would do for the upper bound. It is somewhat interesting to me that the +1SD and -1SD hasn't failed me yet (NCoupons from 2-10,000 tested).

edit: To add the latest version. Put in another ifelse() to assist with a possible error from silly users who input 1 for NCoupons.

Code:

coupon <- function(NCoupons, MEDIAN=TRUE, tol=.001) {
  p_i <- (NCoupons - (seq(1,NCoupons) - 1)) / NCoupons
  EV <- sum(1/p_i)
  VAR <- sum((1-p_i)/p_i^2)
  SD <- sqrt(VAR)  
  if(MEDIAN==TRUE) {
    med.fun <- function(N, NCoupons) {
      seq.max <- ifelse(NCoupons > 100, 100, NCoupons)
      s1 <- seq(1,seq.max)
      x <- sum((-1)^(s1+1) * choose(NCoupons,s1)*((NCoupons-s1)/NCoupons)^N)
      x - .5
    }
    MED <- ifelse(NCoupons==1, 1, uniroot(med.fun, NCoupons=NCoupons, c(EV-SD, EV+SD), tol=tol)$root)
    res <- cbind(EV, MED, VAR, SD)
    rownames(res) <- c("Results")
    colnames(res) <- c("Mean", "MED", "VAR", "SD")
    return(res)  
  }
  res <- cbind(EV, VAR, SD)
  rownames(res) <- c("Results")
  colnames(res) <- c("Mean", "VAR", "SD")
  return(res)  
}

Last edited by Sherman; 12-20-2012 at 02:34 AM.

Quote

12-20-2012 , 02:33 AM

#40

BruceZ

Carpal \'Tunnel

Join Date: Sep 2002 Posts: 11,877

It should be due to non-independence. Nevermind about 23820 failing as a lower bound. The right answer is below that, duh. ds was giving too high a lower bound. I made it lower. This works.

Code:

ds = function(C) {
    lower = ceiling(log(C)/log(2)-1) * floor(0.7*(C-1))
    upper = ceiling(log(C)/log(2)+1) * floor(0.7*(C-1))
    c(lower,upper)  
}

For some reason 195 is now failing with a lower bound of 195 like yours was, and that was working before.

Quote

12-20-2012 , 02:41 AM

#41

Sherman

Carpal \'Tunnel

Join Date: Jun 2005 Posts: 7,766

Quote:

Originally Posted by BruceZ

It should be due to non-independence. Nevermind about 23820 failing as a lower bound. The right answer is below that, duh. ds was giving too high a lower bound. I made it lower. This works.

Code:

ds = function(C) {
    lower = ceiling(log(C)/log(2)-1) * floor(0.7*(C-1))
    upper = ceiling(log(C)/log(2)+1) * floor(0.7*(C-1))
    c(lower,upper)  
}

For some reason 195 is now failing with a lower bound of 195 like yours was, and that was working before.

Now I am getting the version to work that uses ds to get the bounds. It works for 2838 and for 10,000. However, it fails for numbers below 9 (2-8). It is also failing for 15 and 16. Works for all other numbers up to 10,000. Why would 2-8, 15, and 16 not work?

Quote

12-20-2012 , 02:45 AM

#42

BruceZ

Carpal \'Tunnel

Join Date: Sep 2002 Posts: 11,877

It fails below 100 because we're adding 100 terms. You can't add more terms than coupons.

Quote

12-20-2012 , 02:48 AM

#43

Sherman

Carpal \'Tunnel

Join Date: Jun 2005 Posts: 7,766

Quote:

Originally Posted by BruceZ

It fails below 100 because we're adding 100 terms. You can't add more terms than coupons.

Nah. I fixed that problem by adding an ifelse statement. If NCoupons is below 100, then it uses NCoupons as the maximum number of terms.

And it works for every number except 2-8, 15, and 16. Code below for clarity:

Code:

coupon.ds <- function(NCoupons, MEDIAN=TRUE, tol=.001) {
  p_i <- (NCoupons - (seq(1,NCoupons) - 1)) / NCoupons
  EV <- sum(1/p_i)
  VAR <- sum((1-p_i)/p_i^2)
  SD <- sqrt(VAR)  
  if(MEDIAN==TRUE) {
    med.fun <- function(N, NCoupons) {
      seq.max <- ifelse(NCoupons > 100, 100, NCoupons)
      s1 <- seq(1,seq.max)
      x <- sum((-1)^(s1+1) * choose(NCoupons,s1)*((NCoupons-s1)/NCoupons)^N)
      x - .5
    }
    lower <- ceiling(log(NCoupons)/log(2)-1) * floor(0.7*(NCoupons-1))
    upper <- ceiling(log(NCoupons)/log(2)+1) * floor(0.7*(NCoupons-1))
    MED <- uniroot(med.fun, NCoupons=NCoupons, c(lower, upper), tol=tol)$root
    res <- cbind(EV, MED, VAR, SD)
    rownames(res) <- c("Results")
    colnames(res) <- c("Mean", "MED", "VAR", "SD")
    return(res)  
  }
  res <- cbind(EV, VAR, SD)
  rownames(res) <- c("Results")
  colnames(res) <- c("Mean", "VAR", "SD")
  return(res)  
}

Quote

12-20-2012 , 02:51 AM

#44

BruceZ

Carpal \'Tunnel

Join Date: Sep 2002 Posts: 11,877

Replace 100 with min(100,NCoupons). But it still fails at 9 because ds upper bound is too low. It's 16 and needs to be at least 23.

Quote

12-20-2012 , 02:54 AM

#45

BruceZ

Carpal \'Tunnel

Join Date: Sep 2002 Posts: 11,877

Here's a new ds that works all the way down to 1. I just added 10 to the upper bound.

Code:

ds = function(C) {
    lower = ceiling(log(C)/log(2)-1) * floor(0.7*(C-1))
    upper = ceiling(log(C)/log(2)+1) * floor(0.7*(C-1)) + 10
    c(lower,upper)  
}

Quote

12-20-2012 , 02:55 AM

#46

Sherman

Carpal \'Tunnel

Join Date: Jun 2005 Posts: 7,766

Quote:

Originally Posted by BruceZ

Replace 100 with min(100,NCoupons). But it still fails at 9 because ds upper bound is too low. It's 16 and needs to be at least 23.

That is a better way of writing it, but the ifelse() statement I was using did the same thing. If NCoupons > 100, it used 100, else it used NCoupons. So if NCoupons was 100 or less, it used whatever NCoupons was as its number. In other words, the min of 100 or NCoupons.

9 Works for me. It is 8 and fewer that don't. But I see what you mean about why they don't work.

Quote

12-20-2012 , 02:58 AM

#47

Sherman

Carpal \'Tunnel

Join Date: Jun 2005 Posts: 7,766

Quote:

Originally Posted by BruceZ

Here's a new ds that works all the way down to 1. I just added 10 to the upper bound.

Code:

ds = function(C) {
    lower = ceiling(log(C)/log(2)-1) * floor(0.7*(C-1))
    upper = ceiling(log(C)/log(2)+1) * floor(0.7*(C-1)) + 10
    c(lower,upper)  
}

Bingo. That does it. No errors from 1-10,000.

Quote

12-20-2012 , 03:00 AM

#48

Sherman

Carpal \'Tunnel

Join Date: Jun 2005 Posts: 7,766

Side note. This was fun and I learned a lot. Thanks BruceZ!

Quote

12-20-2012 , 03:01 AM

#49

BruceZ

Carpal \'Tunnel

Join Date: Sep 2002 Posts: 11,877

We can do other quantiles next. That should be easy.

Quote

12-20-2012 , 03:01 AM

#50

Sherman

Carpal \'Tunnel

Join Date: Jun 2005 Posts: 7,766

Quote:

Originally Posted by BruceZ

We can do other quantiles next. That should be easy.

Ah yes. I forgot about that suggestion you made earlier. But it is 2am and I am due for bed.

Here are the 2 versions that I currently have:

Code:

coupon <- function(NCoupons, MEDIAN=TRUE, tol=.001) {
  p_i <- (NCoupons - (seq(1,NCoupons) - 1)) / NCoupons
  EV <- sum(1/p_i)
  VAR <- sum((1-p_i)/p_i^2)
  SD <- sqrt(VAR)  
  if(MEDIAN==TRUE) {
    med.fun <- function(N, NCoupons) {
      seq.max <- min(100, NCoupons)
      s1 <- seq(1,seq.max)
      x <- sum((-1)^(s1+1) * choose(NCoupons,s1)*((NCoupons-s1)/NCoupons)^N)
      x - .5
    }
    MED <- ifelse(NCoupons==1, 1, uniroot(med.fun, NCoupons=NCoupons, c(EV-SD, EV+SD), tol=tol)$root)
    res <- cbind(EV, MED, VAR, SD)
    rownames(res) <- c("Results")
    colnames(res) <- c("Mean", "MED", "VAR", "SD")
    return(res)  
  }
  res <- cbind(EV, VAR, SD)
  rownames(res) <- c("Results")
  colnames(res) <- c("Mean", "VAR", "SD")
  return(res)  
}

Code:

coupon.ds <- function(NCoupons, MEDIAN=TRUE, tol=.001) {
  p_i <- (NCoupons - (seq(1,NCoupons) - 1)) / NCoupons
  EV <- sum(1/p_i)
  VAR <- sum((1-p_i)/p_i^2)
  SD <- sqrt(VAR)  
  if(MEDIAN==TRUE) {
    med.fun <- function(N, NCoupons) {
      seq.max <- min(100, NCoupons)
      s1 <- seq(1,seq.max)
      x <- sum((-1)^(s1+1) * choose(NCoupons,s1)*((NCoupons-s1)/NCoupons)^N)
      x - .5
    }
    lower <- ceiling(log(NCoupons)/log(2)-1) * floor(0.7*(NCoupons-1))
    upper <- ceiling(log(NCoupons)/log(2)+1) * floor(0.7*(NCoupons-1)) + 10
    MED <- ifelse(NCoupons==1, 1, uniroot(med.fun, NCoupons=NCoupons, c(lower, upper), tol=tol)$root)
    res <- cbind(EV, MED, VAR, SD)
    rownames(res) <- c("Results")
    colnames(res) <- c("Mean", "MED", "VAR", "SD")
    return(res)  
  }
  res <- cbind(EV, VAR, SD)
  rownames(res) <- c("Results")
  colnames(res) <- c("Mean", "VAR", "SD")
  return(res)  
}

Quote

Page 2 of 4

First

1 2 3 4

Last

Post Reply Subscribe

...

Page 2 of 4

First

1 2 3 4

Last