rggplot2plotbinomial-cdf

Find the two points on the x-axis corresponding to the intersections for a cdf


The corresponding r code is given below.

theta <- seq(0,1, length = 10)
CD_theta <- function(x, p, n){
  1 - pbinom(x, size = n, prob = p) + 1 / 2 * dbinom(x, size = n, prob = p)
}

Then I graphed the data as follows :

mytheta <- CD_theta(5, theta, 10)
df <- data.frame(theta = theta, mytheta = mytheta)
ggplot(df, aes(x = theta, y = mytheta)) +
  geom_line(size = 1, col = "steelblue") +
  ylab("H(Theta)") +
  xlab("Theta")

The corresponding graph is given below. enter image description here

As you can see there are two horizontal lines (graphed in red) and the two vertical lines (graphed in black). I need to find the two points on the x-axis corresponding to the intersections for a H(theta).

I used the locator() function in r to compute the two x intercepts for a single iteration. I would like iterate the above for 1000 times, so it is really tedious approach find them separately.

are there any other r functions can be used find these two x intercept points?

Thank you in advance.


Solution

  • Here is a numerical approach using optimize function:

    library(reprex)
    
    theta <- seq(0,1, length = 10)
    CD_theta <- function(x, p, n){
      1 - pbinom(x, size = n, prob = p) + 1 / 2 * dbinom(x, size = n, prob = p)
    }
    
    # Create a function to compute the difference between the "y" you want 
    # and the "y" computed with CD_theta function 
    # then try to minimize the output of this new function : 
    # get the theta value corresponding to this "y"
    
    my_fn <- function(theta_loc, y, x, n) { 
      # the function to optimize
      abs(y - CD_theta(x, theta_loc, n)) # dont forget abs (absolute)
    }
    
    # Then use optimize function to compute the theta value 
    # between a given interval : c(0,1) in this case
    # Note that you can directly modify here the values of y, x and n
    v1 <- optimize(my_fn, c(0, 1), y = 0.025, x = 5, n = 10)$`minimum` 
    v2 <- optimize(my_fn, c(0, 1), y = 0.975, x = 5, n = 10)$`minimum` 
    
    # print the results
    v1 # 0.025
    #> [1] 0.2120079
    v2 # 0.975
    #> [1] 0.7879756
    

    Created on 2018-09-21 by the reprex package (v0.2.0).