On Jun 20, 2013, at 2:13 PM, Ian Tickle <[email protected]> wrote: > Douglas, I think you are missing the point that estimation of the parameters > of the proper Bayesian statistical model (i.e. the Wilson prior) in order to > perform the integration in the manner you are suggesting requires knowledge > of the already integrated intensities!
Well, that's true, but that's how FW do it. They allow for negative integrated intensities. I'm arguing that we should not do that, since true intensities are positive, and so any estimate of them should also be positive. Examples are always better than words, so here goes (and I apologize for the length): The current way of doing things is summarized by Ed's equation: Ispot-Iback=Iobs. Here Ispot is the # of counts in the spot (the area encompassing the predicted reflection), and Iback is # of counts in the background (usu. some area around the spot). Our job is to estimate the true intensity Itrue. Ed and others argue that Iobs is a reasonable estimate of Itrue, but I say it isn't because Itrue can never be negative, whereas Iobs can. Now where does the Ispot-Iback=Iobs equation come from? It implicitly assumes that both Iobs and Iback come from a Gaussian distribution, in which Iobs and Iback can have negative values. Here's the implicit data model: Ispot = Iobs + Iback There is an Itrue, to which we add some Gaussian noise and randomly generate an Iobs. To that is added some background noise, Iback, which is also randomly generated from a Gaussian with a "true" mean of Ibtrue. This gives us the Ispot, the measured intensity in our spot. Given this data model, Ispot will also have a Gaussian distribution, with mean equal to the sum of Itrue + Ibtrue. From the properties of Gaussians, then, the ML estimate of Itrue will be Ispot-Iback, or Iobs. Now maybe you disagree with that Gaussian data model. If so, welcome to my POV. There are better models, ones that don't give Ispot-Iback as our best estimate of Itrue. Here is a simple example that incorporates our knowledge that Itrue cannot be negative (this example is primarily for illustrating the point, it's not exactly what I would recommend). Instead of using Gaussians, we will use Gamma distributions, which cannot be negative. We assume Iobs is distributed according to a Gamma(Itrue,1). The mean of this distribution is Itrue. (The Maxwell-Boltzmann energy distribution is also a gamma, just for comparison). We also assume that the noise is exponential (a special case of the gamma), Gamma(1,1). The mean of this distribution is 1. (You could imagine that you've normalized Ispot relative to its background --- again, just for ease of calculation). We still assume that Ispot = Iobs + Iback. Then, Ispot will also have a gamma distribution, Gamma(Itrue+1,1). The mean of the Ispot distribution, as you might expect, is Itrue+1. Now we measure Ispot. Given Ispot, the ML estimate of Itrue is: InvDiGamma[ln(Ispot)]-1 if Ispot>0.561 or 0 if Ispot<0.561 Note, the ML estimate is no longer Iobs, and the ML estimate cannot be negative. InvDiGamma is the inverse Digamma function --- a bit unusual, but easily calculated (actually no weirder than the exponential or logarithm, its a relative of factorial and the gamma function). Not something the Braggs would've used, but hey, we've got iPhones now. We can also estimate the SD of of our estimate, but I won't bore you with the equation. A few examples: Ispot ML Itrue SD ----- -------- -- 0.5 0 0.78 0.6 0.04 0.80 0.8 0.25 0.91 0.9 0.36 0.97 1.0 0.46 1.0 1.5 0.97 1.2 2.0 1.48 1.4 3.0 2.49 1.7 5.0 4.49 2.2 10.0 9.50 3.2 20.0 19.5 4.5 100 99.5 10 Note that the first four entries in the table are the case when Ispot<Iback. No negative estimates. You'd get qualitatively similar results if you assume Poisson for Iback and Ispot. To sum up --- the equation Ispot-Iback=Iobs is unphysical because it is founded on unphysical assumptions. If you make better physical assumptions (i.e., Itrue cannot be negative), you end up with different estimates for Itrue. > I suppose we could iterate, i.e. assume an approximate prior, integrate, > calculate a better prior, re-do the integration with the new prior and so on > (hoping of course that the whole process converges), but I think most people > would regard that as overkill. Also dealing with the issue of averaging > estimates of intensities that no longer have a Gaussian error distribution, > and also crucially outlier rejection, would require some rethinking of the > algorithms. The question is would it make any difference in the end compared > with the 'post-correction' we're doing now? > > Cheers > > -- Ian > > > On 20 June 2013 18:14, Douglas Theobald <[email protected]> wrote: > I still don't see how you get a negative intensity from that. It seems you > are saying that in many cases of a low intensity reflection, the integrated > spot will be lower than the background. That is not equivalent to having a > negative measurement (as the measurement is actually positive, and sometimes > things are randomly less positive than backgroiund). If you are using a > proper statistical model, after background correction you will end up with a > positive (or 0) value for the integrated intensity. > > > On Jun 20, 2013, at 1:08 PM, Andrew Leslie <[email protected]> wrote: > > > > > The integration programs report a negative intensity simply because that is > > the observation. > > > > Because of noise in the Xray background, in a large sample of intensity > > estimates for reflections whose true intensity is very very small one will > > inevitably get some measurements that are negative. These must not be > > rejected because this will lead to bias (because some of these intensities > > for symmetry mates will be estimated too large rather than too small). It > > is not unusual for the intensity to remain negative even after averaging > > symmetry mates. > > > > Andrew > > > > > > On 20 Jun 2013, at 11:49, Douglas Theobald <[email protected]> wrote: > > > >> Seems to me that the negative Is should be dealt with early on, in the > >> integration step. Why exactly do integration programs report negative Is > >> to begin with? > >> > >> > >> On Jun 20, 2013, at 12:45 PM, Dom Bellini <[email protected]> > >> wrote: > >> > >>> Wouldnt be possible to take advantage of negative Is to > >>> extrapolate/estimate the decay of scattering background (kind of Wilson > >>> plot of background scattering) to flat out the background and push all > >>> the Is to positive values? > >>> > >>> More of a question rather than a suggestion ... > >>> > >>> D > >>> > >>> > >>> > >>> From: CCP4 bulletin board [mailto:[email protected]] On Behalf Of Ian > >>> Tickle > >>> Sent: 20 June 2013 17:34 > >>> To: ccp4bb > >>> Subject: Re: [ccp4bb] ctruncate bug? > >>> > >>> Yes higher R factors is the usual reason people don't like I-based > >>> refinement! > >>> > >>> Anyway, refining against Is doesn't solve the problem, it only postpones > >>> it: you still need the Fs for maps! (though errors in Fs may be less > >>> critical then). > >>> -- Ian > >>> > >>> On 20 June 2013 17:20, Dale Tronrud > >>> <[email protected]<mailto:[email protected]>> wrote: > >>> If you are refining against F's you have to find some way to avoid > >>> calculating the square root of a negative number. That is why people > >>> have historically rejected negative I's and why Truncate and cTruncate > >>> were invented. > >>> > >>> When refining against I, the calculation of (Iobs - Icalc)^2 couldn't > >>> care less if Iobs happens to be negative. > >>> > >>> As for why people still refine against F... When I was distributing > >>> a refinement package it could refine against I but no one wanted to do > >>> that. The "R values" ended up higher, but they were looking at R > >>> values calculated from F's. Of course the F based R values are lower > >>> when you refine against F's, that means nothing. > >>> > >>> If we could get the PDB to report both the F and I based R values > >>> for all models maybe we could get a start toward moving to intensity > >>> refinement. > >>> > >>> Dale Tronrud > >>> > >>> > >>> On 06/20/2013 09:06 AM, Douglas Theobald wrote: > >>> Just trying to understand the basic issues here. How could refining > >>> directly against intensities solve the fundamental problem of negative > >>> intensity values? > >>> > >>> > >>> On Jun 20, 2013, at 11:34 AM, Bernhard Rupp > >>> <[email protected]<mailto:[email protected]>> wrote: > >>> As a maybe better alternative, we should (once again) consider to refine > >>> against intensities (and I guess George Sheldrick would agree here). > >>> > >>> I have a simple question - what exactly, short of some sort of historic > >>> inertia (or memory lapse), is the reason NOT to refine against > >>> intensities? > >>> > >>> Best, BR > >>> > >>> > >>> > >>> > >>> -- > >>> > >>> This e-mail and any attachments may contain confidential, copyright and > >>> or privileged material, and are for the use of the intended addressee > >>> only. If you are not the intended addressee or an authorised recipient of > >>> the addressee please notify us of receipt by returning the e-mail and do > >>> not use, copy, retain, distribute or disclose the information in or > >>> attached to the e-mail. > >>> > >>> Any opinions expressed within this e-mail are those of the individual and > >>> not necessarily of Diamond Light Source Ltd. > >>> > >>> Diamond Light Source Ltd. cannot guarantee that this e-mail or any > >>> attachments are free from viruses and we cannot accept liability for any > >>> damage which you may sustain as a result of software viruses which may be > >>> transmitted in or with the message. > >>> > >>> Diamond Light Source Limited (company no. 4375679). Registered in England > >>> and Wales with its registered office at Diamond House, Harwell Science > >>> and Innovation Campus, Didcot, Oxfordshire, OX11 0DE, United Kingdom > >>> > >>> > >>> > >>> > >>> > >>> > >>> > >>> > >>> > > >
