title: Simulation for Confidence interval in R
body: <p>I have an R function that provides the 95% confidence Interval for the <em>ncp</em> (non-centrality parameter) of a t distribution.</p>
  
  <p>Via simulation in R, is it possible to show that in the long-run the CIs from this R function capture a given <em>TRUE ncp</em> (here "2" same as input <em>t</em>) 95%  of the time? </p>
  
  <p>(I appreciate any ideas as to how to do this)</p>
  
  <pre><code>CI.ncp &lt;- function(t, N){
  
    f &lt;- function (ncp, alpha, q, df) {  
  abs(suppressWarnings(pt(q = t, df = N - 1, ncp, lower.tail = FALSE)) - alpha) }
  
  sapply(c(0.025, 0.975),
  function(x) optim(1, f, alpha = x, q = t, df = N - 1, control = list(reltol = (.Machine$double.eps)))[[1]]) }
  
  #Example of Use:
  CI.ncp(t = 2, N = 20) # gives: -0.08293755  4.03548862 
  
  #(in the long-run 95% of the time, "2" is contained within these
  # two numbers, how to show this in R?)
  </code></pre>
  
  <p>Here is what I have tried with no success:</p>
  
  <pre><code>fun &lt;- function(t = 2, N = 20){
  
    ncp = rt(1, N - 1, t)
    CI.ncp(t = 2, N = 20)
    mean(ncp &lt;= 2 &amp; 2 &lt;= ncp )
     }
  
   R &lt;- 1000
   sim &lt;- t(replicate(R, fun()))
   coverage &lt;- mean(sim[,1] &lt;= 2 &amp; 2 &lt;= sim[,2])
  </code></pre>