4085
/******************* Fast Lambda Calculation Subroutine ************************
4086
Version 1.0 May 16, 1991
4087
Program by: Stephen Altschul
4089
Uses Newton-Raphson method (fast) to solve for Lambda, given an initial
4090
guess (lambda0) obtained perhaps by the bisection method.
4091
*******************************************************************************/
4094
* Find positive solution to sum_{i=low}^{high} exp(i lambda) = 1.
4096
* @param probs probabilities of a score occurring
4097
* @param d the gcd of the possible scores. This equals 1 if the scores
4099
* @param low the lowest possible score
4100
* @param high the highest possible score
4101
* @param lambda0 an initial value for lambda
4102
* @param tolx the tolerance to which lambda must be computed
4103
* @param itmax the maximum number of times the function may be
4105
* @param maxNewton the maximum permissible number of Newton
4106
* iteration. After that the computation will proceed by bisection.
4107
* @param itn a pointer to an integer that will receive the actually
4108
* number of iterations performed.
4110
* Let phi(lambda) = sum_{i=low}^{high} exp(i lambda) - 1. Then phi(lambda)
4113
* phi(lamdba) = exp(u lambda) p( exp(-lambda) )
4115
* where p(x) is a polynomial that has exactly two zeros, one at x = 1
4116
* and one at y = exp(-lamdba). It is simpler, more numerically
4117
* efficient and stable to apply Newton's method to p(x) than to
4120
* We define a safeguarded Newton iteration as follows. Let the
4121
* initial interval of uncertainty be [0,1]. If p'(x) >= 0, we bisect
4122
* the interval. Otherwise we try a Newton step. If the Newton iterate
4123
* lies in the current interval of uncertainty and it reduces the
4124
* value of | p(x) | by at least 10%, we accept the new
4125
* point. Otherwise, we bisect the current interval of uncertainty.
4126
* It is clear that this method converges to a zero of p(x). Since
4127
* p'(x) > 0 in an interval containing x = 1, the method cannot
4128
* converge to x = 1 and therefore converges to the only other zero,
4089
* Find positive solution to
4091
* sum_{i=low}^{high} exp(i lambda) * probs[i] = 1.
4093
* Note that this solution does not exist unless the average score is
4094
* negative and the largest score that occurs with nonzero probability
4097
* @param probs probabilities of a score occurring
4098
* @param d the gcd of the possible scores. This equals 1 if
4099
* the scores are not a lattice
4100
* @param low the lowest possible score that occurs with
4101
* nonzero probability
4102
* @param high the highest possible score that occurs with
4103
* nonzero probability.
4104
* @param lambda0 an initial guess for lambda
4105
* @param tolx the tolerance to which lambda must be computed
4106
* @param itmax the maximum number of times the function may be
4108
* @param maxNewton the maximum permissible number of Newton
4109
* iterations; after that the computation will proceed
4111
* @param *itn the number of iterations needed to compute Lambda,
4112
* or itmax if Lambda could not be computed.
4114
* Let phi(lambda) = sum_{i=low}^{high} exp(i lambda) - 1. Then
4115
* phi(lambda) may be written
4117
* phi(lamdba) = exp(u lambda) f( exp(-lambda) )
4119
* where f(x) is a polynomial that has exactly two zeros, one at x = 1
4120
* and one at x = exp(-lamdba). It is simpler to solve this problem
4121
* in x = exp(-lambda) than it is to solve it in lambda, because we
4122
* know that for x, a solution lies in [0,1], and because Newton's
4123
* method is generally more stable and efficient for polynomials than
4124
* it is for exponentials.
4126
* For the most part, this function is a standard safeguarded Newton
4127
* iteration: define an interval of uncertainty [a,b] with f(a) > 0
4128
* and f(b) < 0 (except for the initial value b = 1, where f(b) = 0);
4129
* evaluate the function and use the sign of that value to shrink the
4130
* interval of uncertainty; compute a Newton step; and if the Newton
4131
* step suggests a point outside the interval of uncertainty or fails
4132
* to decrease the function sufficiently, then bisect. There are
4133
* three further details needed to understand the algorithm:
4135
* 1) If y the unique solution in [0,1], then f is positive to the left of
4136
* y, and negative to the right. Therefore, we may determine whether
4137
* the Newton step -f(x)/f'(x) is moving toward, or away from, y by
4138
* examining the sign of f'(x). If f'(x) >= 0, we bisect instead
4139
* of taking the Newton step.
4140
* 2) There is a neighborhood around x = 1 for which f'(x) >= 0, so
4141
* (1) prevents convergence to x = 1 (and for a similar reason
4142
* prevents convergence to x = 0, if the function is incorrectly
4143
* called with probs[high] == 0).
4144
* 3) Conditions like fabs(p) < lambda_tolerance * x * (1-x) are used in
4145
* convergence criteria because these values translate to a bound
4146
* on the relative error in lambda. This is proved in the
4147
* "Blast Scoring Parameters" document that accompanies the BLAST
4150
* The iteration on f(x) is robust and doesn't overflow; defining a
4151
* robust safeguarded Newton iteration on phi(lambda) that cannot
4152
* converge to lambda = 0 and that is protected against overflow is
4153
* more difficult. So (despite the length of this comment) the Newton
4154
* iteration on f(x) is the simpler solution.
4132
4157
static Nlm_FloatHi
4133
NlmKarlinLambdaNR( Nlm_FloatHi PNTR probs, BLAST_Score d,
4134
BLAST_Score low, BLAST_Score high,
4135
Nlm_FloatHi lambda0, Nlm_FloatHi tolx,
4136
int itmax, int maxNewton, int * itn )
4158
NlmKarlinLambdaNR(Nlm_FloatHi PNTR probs, BLAST_Score d, BLAST_Score low,
4159
BLAST_Score high, Nlm_FloatHi lambda0, Nlm_FloatHi tolx,
4160
int itmax, int maxNewton, int * itn)
4139
4163
Nlm_FloatHi x0, x, a = 0, b = 1;
added
removed
tools/blastkar.c
