6

I'm simulating an extremely rare event (the detection of weighted photon packets in a highly absorbing material). For example, I may simulate the transmission of 1e9 photons and only detect 10 of them (and their weights vary from $0^+$ to 1). I'm calculating the mean value of this occurring (what is the mean of $\sum w_{received}/N_{transmitted}$, where $w_{received}$ is the weight of the received photons and $N_{transmitted}$ is the number of transmitted photons) and put a confidence interval on this. My understanding is that the population mean, $\mu$ lies within one standard Normal deviation of $\bar{x} \pm S_N/\sqrt{N}$, where $\bar{x}$ is the sample mean, $S_N$ is the sample variance, and $N$ is the number of samples.

My question is this, how do I know when I've received enough "positive" samples, or received photons, to get an accurate estimate? Since I'm dealing with such a small probability, the total number of samples, $N$ is extremely high and $S_N/\sqrt{N}$ will always be small. I've seen notes saying the CLT only holds for samples sizes > 30, but would that apply to something like this? i.e. $N_{received}$ > 30?

Summary: I'm drawing samples from an unknown distribution and computing their mean. Sample values range from 0 to 1, with 0 being far more common. I'm computing the mean (which is very close to 0) and I want to put a confidence interval on this.

Edit From my description above, there's two ways to think about the problem. The first is a binomial distribution of yes/no results (yes, received; no, not received). However, the more complex, and in practice, more useful, distribution is somewhat of a weighted binomial (not sure if that's the right term), where each received photon "packet" has a power associated with it, that ranges from 0 to 1. That is the case I care most about putting a confidence interval on. i.e. my sample values will predominately be 0's with a small number of samples with values 0 < x <= 1

Edit 2 The PDF will look something like the crude drawing below. Note that it is not to scale - there will be much more 0-weights than non-zero weights. This is an example pdf of samples. I don't know what this will look like before I do an experiment - I just know my samples will range from 0 to 1, with a predominate amount being 0.

Example PDF of sample values

Edit 3 The term "weight" seems to be a bit misleading. "Sample value" or "power" might be a better term. My simulation tracks photon movement through water. Most photons never reach the receiver and therefore have a power of 0. As photons move through the water their power is reduced (to model a group of photons moving with a portion being absorbed as they move along). As these photons are received they have weights that range from a max of 1 (no loss from the travel) down towards a number approaching 0.

gallamine
  • 373
  • 1
    The mean value of the observed frequencies is an inferior estimator when $N_{\text{transmitted}}$ varies appreciably, because the frequencies for the smaller $N$ are less reliable. What matters (usually) is the total received divided by the total transmitted. Poisson confidence limits will work well. – whuber Sep 15 '11 at 05:48
  • Thanks whuber, $N_{transmitted}$ is usually fixed and $N_{received}$ varies. – gallamine Sep 15 '11 at 12:31
  • Re the edit: What exactly does the "weight" mean? Does it mean that the probabilities of detection are varying from one packet to the next? If so, are these probabilities changing in known ways or unknown ways? – whuber Sep 15 '11 at 16:16
  • Whuber, sorry I guess it can be a bit difficult to explain. By "weight" I mean the value of the detected event/photon. No detection yields a weight of 0 and detection yields a weight from $0^+$ to 1. So, I launch 1E9 photons, most are terminated (weight = 0) while a small portion are detected with varying weights. I'm trying to find the mean value of weights along with a confidence interval. – gallamine Sep 15 '11 at 17:57
  • What is a detection weight? Is it perhaps an uncertain detection, or is it something else? Note that your problem has now fundamentally changed: it is no longer about a Binomial (or Poisson) distribution and requires a different approach. – whuber Sep 15 '11 at 18:44
  • 1
    Whuber, Yes I recognize that it changes things and I should have been more explicit in my description.

    To speed computation time, it's customary to treat ray tracing problems of photon movement or neutron movement as groups of photons instead of individual photons. That way, instead of terminating the trajectory at an absorption event, the group weight is simply reduced. The computation time is not lost and the aggregate received power is measured. A weight of 1 indicates no photons in the "group" was lost due to absorption. Weight=0 is non-detected and detected photons have weights=(0,1]

     – gallamine Sep 15 '11 at 20:13

2 Answers2

4

The normal approximation for the confidence interval of binomial proportions breaks down very badly for rare events and the rules of thumb about sample sizes are inconsistent and unreliable. Better methods are just as easy to calculate (i.e. you click the button!) and so there is no reason for anyone to use the normal approximation. Ever.

Have a quick look at the papers below (and then use Wilson's method).

Vollset. Confidence intervals for a binomial proportion. Statist. Med. (1993) vol. 12 (9) pp. 809-24 Brown et al. Interval Estimation for a Binomial Proportion. Statistical Science (2001) pp. 101-117 http://www.jstor.org/stable/2676784

See also some previous questions put to this list: How to report asymmetrical confidence intervals of a proportion? and Discrete functions: Confidence interval coverage? and Clarification on interpreting confidence intervals?

Community
  • 1
Michael Lew
  • 15,102
  • I believe a previous answer, which is a good one, sharply disagrees with your prefatory comments: http://stats.stackexchange.com/questions/4756/confidence-interval-for-bernoulli-sampling/6184#6184. Wilson's method looks like a good recommendation, though. – whuber Sep 15 '11 at 05:53
  • Michael, thanks for the answer. I'm reading through your provided material. However, I do wonder if the answer applies when the sample distribution isn't strictly binomial, but weighted? I edited my question to provide more details. Essentially, in my case, I either don't receive a photon packet (weight = 0) or I do, and the weight varies (0 < weight <= 1). – gallamine Sep 15 '11 at 12:28
  • 1
    Gallamine, I don't really know anything about intervals for weighted binomial distributions. However, I would be surprised if a weighting made the normal approximation any more reliable. – Michael Lew Sep 16 '11 at 04:03
  • 1
    Whuber, I am happy to be in disagreement with the answer that you point to: the papers that I cited are just some of many that show erratic coverage failures for intervals from the normal approximation, and better intervals are not difficult to determine. – Michael Lew Sep 16 '11 at 04:07
1

Now that it is clear that you have a weighting function, I suggest that you use Bayesian intervals (often called credible intervals) with the weighting function being the prior. Multiply that by the likelihood function provided by your results to get the posterior. Any interval containing 95% of the area under that posterior distribution is a 95% credible interval.

The likelihood function is easily calculated: start with a uniform (0,1) indicating no data and so no evidence. For each photon received you multiply the distribution by y=x and for each photon sent but not received multiply it by y=1-x. When you've done that for all of the photons sent you will have the likelihood function representing the evidence inherent in your data. You can scale it to a maximum of 1 to look conventional, if you like. [Of course, y represents the likelihood and x is the hypothetical probability of success in each trial.]

There is a formula for the likelihood function, but I find it easier to understand in the way I've expressed it here.

Michael Lew
  • 15,102
  • Michael, I confess that I'm not quite certain what you're getting at here. My "weighting function", or the distribution of sample values is unknown prior to my simulation. I only know this after the fact. The image in my question is just an example of what it might look like. In general, I just know that my values will be between 0 and 1. – gallamine Sep 19 '11 at 16:07
  • Gallamine, it seems possible to me that when you say PDF you really mean likelihood function... Conventionally the confidence interval is for the true (fixed) probability of success for each trial. The notion of a PDF in that circumstance fit well with the idea that you had information about the distribution of possible values for that true probability of success... Oh well, back to the Wilson's interval I guess. – Michael Lew Sep 20 '11 at 04:06