1- using Numerics . Mathematics . Optimization ;
1+ using Numerics . Mathematics . Optimization ;
22using System ;
33using System . Collections . Generic ;
44
@@ -24,17 +24,35 @@ public class YeoJohnson
2424 /// <param name="lambda">Output. The transformation exponent. Range -5 to +5.</param>
2525 public static void FitLambda ( IList < double > values , out double lambda )
2626 {
27- // Solve with Brent
28- var brent = new BrentSearch ( ( x ) => { return LogLikelihood ( values , x ) ; } , - 5d , 5d ) ;
29- brent . Maximize ( ) ;
30- if ( brent . Status == OptimizationStatus . Success )
27+ lambda = double . NaN ;
28+ if ( ! CanFitLambda ( values ) )
29+ return ;
30+
31+ // Keep BrentSearch arithmetic finite even when the profile likelihood is undefined.
32+ var brent = new BrentSearch ( ( x ) =>
3133 {
32- lambda = brent . BestParameterSet . Values [ 0 ] ;
33- }
34- else
34+ double logLikelihood = LogLikelihood ( values , x ) ;
35+ return Tools . IsFinite ( logLikelihood ) ? logLikelihood : - double . MaxValue ;
36+ } , - 5d , 5d )
3537 {
36- lambda = double . NaN ;
37- }
38+ ReportFailure = false ,
39+ ComputeHessian = false ,
40+ RecordTraces = false
41+ } ;
42+
43+ brent . Maximize ( ) ;
44+ if ( brent . Status != OptimizationStatus . Success ||
45+ brent . BestParameterSet . Values == null ||
46+ brent . BestParameterSet . Values . Length == 0 )
47+ return ;
48+
49+ double candidate = brent . BestParameterSet . Values [ 0 ] ;
50+ if ( ! Tools . IsFinite ( candidate ) ||
51+ Math . Abs ( candidate ) > 5d ||
52+ ! Tools . IsFinite ( LogLikelihood ( values , candidate ) ) )
53+ return ;
54+
55+ lambda = candidate ;
3856 }
3957
4058 /// <summary>
@@ -62,6 +80,9 @@ public static double FitLambda(IList<double> values)
6280 /// </returns>
6381 public static double LogLikelihood ( IList < double > values , double lambda )
6482 {
83+ if ( ! CanFitLambda ( values ) || ! Tools . IsFinite ( lambda ) || Math . Abs ( lambda ) > 5d )
84+ return double . NegativeInfinity ;
85+
6586 int n = values . Count ;
6687 var transformed = new double [ n ] ;
6788 double sum = 0d ;
@@ -72,8 +93,13 @@ public static double LogLikelihood(IList<double> values, double lambda)
7293 {
7394 double xi = values [ i ] ;
7495 double yi = Transform ( xi , lambda ) ;
96+ if ( ! Tools . IsFinite ( yi ) )
97+ return double . NegativeInfinity ;
98+
7599 transformed [ i ] = yi ;
76100 sum += yi ;
101+ if ( ! Tools . IsFinite ( sum ) )
102+ return double . NegativeInfinity ;
77103
78104 // Compute derivative dT/dy for log-Jacobian
79105 double dTdy ;
@@ -91,6 +117,8 @@ public static double LogLikelihood(IList<double> values, double lambda)
91117 logJacobianSum += Math . Log ( dTdy ) ;
92118 else
93119 return double . NegativeInfinity ; // log-likelihood undefined
120+ if ( ! Tools . IsFinite ( logJacobianSum ) )
121+ return double . NegativeInfinity ;
94122 }
95123
96124 // Compute mean and SSE
@@ -101,14 +129,16 @@ public static double LogLikelihood(IList<double> values, double lambda)
101129 double resid = transformed [ i ] - mu ;
102130 sse += resid * resid ;
103131 }
132+ if ( ! Tools . IsFinite ( sse ) || sse <= 0d )
133+ return double . NegativeInfinity ;
104134
105135 double sigmaSq = sse / n ;
106- if ( sigmaSq <= 0 || double . IsNaN ( sigmaSq ) )
136+ if ( ! Tools . IsFinite ( sigmaSq ) || sigmaSq <= 0 )
107137 return double . NegativeInfinity ;
108138
109139 double logLikelihood = - n / 2.0 * Tools . LogSqrt2PI - n / 2.0 * Math . Log ( sigmaSq ) - sse / ( 2.0 * sigmaSq ) + logJacobianSum ;
110140
111- return logLikelihood ;
141+ return Tools . IsFinite ( logLikelihood ) ? logLikelihood : double . NegativeInfinity ;
112142 }
113143
114144 /// <summary>
@@ -160,11 +190,11 @@ public static double Transform(double value, double lambda)
160190 {
161191 return Math . Log ( value + 1 ) ;
162192 }
163- else if ( value < 0 && lambda != 2 )
193+ else if ( value < 0 && Math . Abs ( lambda - 2d ) >= 1E-8 )
164194 {
165195 return - ( Math . Pow ( - value + 1 , 2 - lambda ) - 1 ) / ( 2 - lambda ) ;
166196 }
167- else if ( value < 0 && lambda == 2 )
197+ else if ( value < 0 && Math . Abs ( lambda - 2d ) < 1E-8 )
168198 {
169199 return - Math . Log ( - value + 1 ) ;
170200 }
@@ -216,11 +246,11 @@ public static double InverseTransform(double value, double lambda)
216246 {
217247 return Math . Exp ( value ) - 1 ;
218248 }
219- else if ( value < 0 && lambda != 2 )
249+ else if ( value < 0 && Math . Abs ( lambda - 2d ) >= 1E-8 )
220250 {
221251 return 1 - ( Math . Pow ( 1 - ( 2 - lambda ) * value , 1 / ( 2 - lambda ) ) ) ;
222252 }
223- else if ( value < 0 && lambda == 2 )
253+ else if ( value < 0 && Math . Abs ( lambda - 2d ) < 1E-8 )
224254 {
225255 return 1 - Math . Exp ( - value ) ;
226256 }
@@ -242,5 +272,32 @@ public static List<double> InverseTransform(IList<double> values, double lambda)
242272 newValues . Add ( InverseTransform ( values [ i ] , lambda ) ) ;
243273 return newValues ;
244274 }
275+
276+ /// <summary>
277+ /// Determines whether a sample can support Yeo-Johnson lambda fitting.
278+ /// </summary>
279+ /// <param name="values">The sample values to inspect.</param>
280+ /// <returns><see langword="true"/> when the sample is finite and non-degenerate.</returns>
281+ private static bool CanFitLambda ( IList < double > values )
282+ {
283+ if ( values == null || values . Count < 2 )
284+ return false ;
285+
286+ double first = values [ 0 ] ;
287+ if ( ! Tools . IsFinite ( first ) )
288+ return false ;
289+
290+ bool hasDifferentValue = false ;
291+ for ( int i = 1 ; i < values . Count ; i ++ )
292+ {
293+ double value = values [ i ] ;
294+ if ( ! Tools . IsFinite ( value ) )
295+ return false ;
296+ if ( value != first )
297+ hasDifferentValue = true ;
298+ }
299+
300+ return hasDifferentValue ;
301+ }
245302 }
246303}
0 commit comments