4
from c_numpy cimport ndarray, npy_intp, \
5
PyArray_SIZE, PyArray_EMPTY, PyArray_FROMANY, \
6
NPY_INT, NPY_DOUBLE, NPY_OWNDATA, NPY_ALIGNED, NPY_FORTRAN, \
7
PyArray_SimpleNewFromData
12
# NumPy must be initialized
13
c_numpy.import_array()
17
cdef floatarray_from_data(int rows, int cols, double *data):
21
a_ndr = <object>PyArray_SimpleNewFromData(1, &size, NPY_DOUBLE, data)
23
a_ndr.shape = (rows, cols)
26
cdef boolarray_from_data(int rows, int cols, int *data):
30
a_ndr = <object>PyArray_SimpleNewFromData(1, &size, NPY_DOUBLE, data)
32
a_ndr.shape = (rows, cols)
33
return a_ndr.astype(numpy.bool)
40
The (m,p) array of independent variables where the surface must be estimated.
42
The (m,) ndarray of loess values evaluated at newdata
44
The (m,) ndarray of the estimates of the standard error on the estimated
46
residual_scale : float
47
Estimate of the scale of the residuals
49
Degrees of freedom of the t-distribution used to compute pointwise
50
confidence intervals for the evaluated surface.
52
Number of new observations.
58
#####---------------------------------------------------------------------------
59
#---- ---- loess model ---
60
#####---------------------------------------------------------------------------
61
cdef class loess_inputs:
66
A (n,p) ndarray of independent variables, with n the number of observations
67
and p the number of variables.
69
A (n,) ndarray of observations
71
Number of observations: nobs=n.
73
Number of independent variables: npar=p.
75
A (n,) ndarray of weights to be given to individual observations in the
76
sum of squared residuals that forms the local fitting criterion. If not
77
None, the weights should be non negative. If the different observations
78
have non-equal variances, the weights should be inversely proportional
80
By default, an unweighted fit is carried out (all the weights are one).
82
cdef c_loess.c_loess_inputs *_base
84
cdef readonly ndarray x, y, masked, x_eff, y_eff
85
cdef readonly int nobs, npar
87
def __init__(self, x_data, y_data):
89
self.x = narray(x_data, copy=False, subok=True, dtype=float_, order='C')
90
self.y = narray(y_data, copy=False, subok=True, dtype=float_, order='C')
91
# Check the dimensions ........
94
elif self.x.ndim == 2:
95
self.npar = self.x.shape[-1]
97
raise ValueError("The array of indepedent varibales should be 2D at most!")
98
self.nobs = len(self.x)
99
# Get the effective data ......
100
self.x_eff = self.x.ravel()
102
self.w_ndr = numpy.ones((self.nobs,), dtype=float_)
106
"""A (n,) ndarray of weights to be given to individual observations in the
107
sum of squared residuals that forms the local fitting criterion. If not
108
None, the weights should be non negative. If the different observations
109
have non-equal variances, the weights should be inversely proportional
111
By default, an unweighted fit is carried out (all the weights are one).
116
def __set__(self, w):
119
w_ndr = narray(w, copy=False, subok=False)
120
if w_ndr.ndim > 1 or w_ndr.size != self.nobs:
121
raise ValueError, "Invalid size of the 'weights' vector!"
123
self._base.weights = <double *>w_ndr.data
125
######---------------------------------------------------------------------------
126
##---- ---- loess control ---
127
######---------------------------------------------------------------------------
128
cdef class loess_control:
129
"""Loess control parameters.
132
surface : string ["interpolate"]
133
Determines whether the fitted surface is computed directly at all points
134
("direct") or whether an interpolation method is used ("interpolate").
135
The default ("interpolate") is what most users should use unless special
136
circumstances warrant.
137
statistics : string ["approximate"]
138
Determines whether the statistical quantities are computed exactly
139
("exact") or approximately ("approximate"). "exact" should only be used
140
for testing the approximation in statistical development and is not meant
141
for routine usage because computation time can be horrendous.
142
trace_hat : string ["wait.to.decide"]
143
Determines how the trace of the hat matrix should be computed. The hat
144
matrix is used in the computation of the statistical quantities.
145
If "exact", an exact computation is done; this could be slow when the
146
number of observations n becomes large. If "wait.to.decide" is selected,
147
then a default is "exact" for n < 500 and "approximate" otherwise.
148
This option is only useful when the fitted surface is interpolated. If
149
surface is "exact", an exact computation is always done for the trace.
150
Setting trace_hat to "approximate" for large dataset will substantially
151
reduce the computation time.
153
Number of iterations of the robust fitting method. If the family is
154
"gaussian", the number of iterations is set to 0.
156
Maximum cell size of the kd-tree. Suppose k = floor(n*cell*span),
157
where n is the number of observations, and span the smoothing parameter.
158
Then, a cell is further divided if the number of observations within it
159
is greater than or equal to k. This option is only used if the surface
163
cdef c_loess.c_loess_control *_base
167
surface : string ["interpolate"]
168
Determines whether the fitted surface is computed directly at all points
169
("direct") or whether an interpolation method is used ("interpolate").
170
The default ("interpolate") is what most users should use unless special
171
circumstances warrant.
174
return self._base.surface
175
def __set__(self, surface):
176
if surface.lower() not in ('interpolate', 'direct'):
177
raise ValueError("Invalid value for the 'surface' argument: "+
178
"should be in ('interpolate', 'direct').")
179
tmpx = surface.lower()
180
self._base.surface = tmpx
184
statistics : string ["approximate"]
185
Determines whether the statistical quantities are computed exactly
186
("exact") or approximately ("approximate"). "exact" should only be used
187
for testing the approximation in statistical development and is not meant
188
for routine usage because computation time can be horrendous.
191
return self._base.statistics
192
def __set__(self, statistics):
193
if statistics.lower() not in ('approximate', 'exact'):
194
raise ValueError("Invalid value for the 'statistics' argument: "\
195
"should be in ('approximate', 'exact').")
196
tmpx = statistics.lower()
197
self._base.statistics = tmpx
201
trace_hat : string ["wait.to.decide"]
202
Determines how the trace of the hat matrix should be computed. The hat
203
matrix is used in the computation of the statistical quantities.
204
If "exact", an exact computation is done; this could be slow when the
205
number of observations n becomes large. If "wait.to.decide" is selected,
206
then a default is "exact" for n < 500 and "approximate" otherwise.
207
This option is only useful when the fitted surface is interpolated. If
208
surface is "exact", an exact computation is always done for the trace.
209
Setting trace_hat to "approximate" for large dataset will substantially
210
reduce the computation time.
213
return self._base.trace_hat
214
def __set__(self, trace_hat):
215
if trace_hat.lower() not in ('approximate', 'exact'):
216
raise ValueError("Invalid value for the 'trace_hat' argument: "\
217
"should be in ('approximate', 'exact').")
218
tmpx = trace_hat.lower()
219
self._base.trace_hat = tmpx
224
Number of iterations of the robust fitting method. If the family is
225
"gaussian", the number of iterations is set to 0.
228
return self._base.iterations
229
def __set__(self, iterations):
231
raise ValueError("Invalid number of iterations: should be positive")
232
self._base.iterations = iterations
237
Maximum cell size of the kd-tree. Suppose k = floor(n*cell*span),
238
where n is the number of observations, and span the smoothing parameter.
239
Then, a cell is further divided if the number of observations within it
240
is greater than or equal to k. This option is only used if the surface
244
return self._base.cell
245
def __set__(self, cell):
247
raise ValueError("Invalid value for the cell argument: should be positive")
248
self._base.cell = cell
250
def update(self, **cellargs):
251
"""Updates several parameters at once."""
252
surface = cellargs.get('surface', None)
253
if surface is not None:
254
self.surface = surface
256
statistics = cellargs.get('statistics', None)
257
if statistics is not None:
258
self.statistics = statistics
260
trace_hat = cellargs.get('trace_hat', None)
261
if trace_hat is not None:
262
self.trace_hat = trace_hat
264
iterations = cellargs.get('iterations', None)
265
if iterations is not None:
266
self.iterations = iterations
268
cell = cellargs.get('cell', None)
270
self.parametric_flags = cell
275
"Surface type : %s" % self.surface,
276
"Statistics : %s" % self.statistics,
277
"Trace estimation : %s" % self.trace_hat,
278
"Cell size : %s" % self.cell,
279
"Nb iterations : %s" % self.iterations,]
280
return '\n'.join(strg)
283
######---------------------------------------------------------------------------
284
##---- ---- loess kd_tree ---
285
######---------------------------------------------------------------------------
286
cdef class loess_kd_tree:
287
cdef c_loess.c_loess_kd_tree *_base
289
######---------------------------------------------------------------------------
290
##---- ---- loess model ---
291
######---------------------------------------------------------------------------
292
cdef class loess_model:
293
"""loess_model contains parameters required for a loess fit.
296
normalize : boolean [True]
297
Determines whether the independent variables should be normalized.
298
If True, the normalization is performed by setting the 10% trimmed
299
standard deviation to one. If False, no normalization is carried out.
300
This option is only useful for more than one variable. For spatial
301
coordinates predictors or variables with a common scale, it should be
304
Smoothing factor, as a fraction of the number of points to take into
307
Overall degree of locally-fitted polynomial. 1 is locally-linear
308
fitting and 2 is locally-quadratic fitting. Degree should be 2 at most.
309
family : string ["gaussian"]
310
Determines the assumed distribution of the errors. The values are
311
"gaussian" or "symmetric". If "gaussian" is selected, the fit is
312
performed with least-squares. If "symmetric" is selected, the fit
313
is performed robustly by redescending M-estimators.
314
parametric_flags : sequence [ [False]*p ]
315
Indicates which independent variables should be conditionally-parametric
316
(if there are two or more independent variables). The argument should be
317
a sequence of booleans, with the same size as the number of independent
318
variables, specified in the order of the predictor group ordered in x.
319
Note: elements of the sequence cannot be modified individually: the whole
320
sequence must be given.
321
drop_square : sequence [ [False]* p]
322
When there are two or more independent variables and when a 2nd order
323
polynomial is used, "drop_square_flags" specifies those numeric predictors
324
whose squares should be dropped from the set of fitting variables.
325
The method of specification is the same as for parametric.
326
Note: elements of the sequence cannot be modified individually: the whole
327
sequence must be given.
331
cdef c_loess.c_loess_model *_base
334
cdef setup(self, c_loess.c_loess_model *base, long npar):
341
return bool(self._base.normalize)
342
def __set__(self, normalize):
343
self._base.normalize = normalize
347
return self._base.span
348
def __set__(self, span):
349
if span <= 0. or span > 1.:
350
raise ValueError("Span should be between 0 and 1!")
351
self._base.span = span
355
return self._base.degree
356
def __set__(self, degree):
357
if degree < 0 or degree > 2:
358
raise ValueError("Degree should be between 0 and 2!")
359
self._base.degree = degree
363
return self._base.family
364
def __set__(self, family):
365
if family.lower() not in ('symmetric', 'gaussian'):
366
raise ValueError("Invalid value for the 'family' argument: "\
367
"should be in ('symmetric', 'gaussian').")
368
self._base.family = family
370
property parametric_flags:
372
return boolarray_from_data(self.npar, 1, self._base.parametric)
373
def __set__(self, paramf):
376
p_ndr = numpy.atleast_1d(narray(paramf, copy=False, subok=True,
378
for i from 0 <= i < min(self.npar, p_ndr.size):
379
self._base.parametric[i] = p_ndr[i]
381
property drop_square_flags:
383
return boolarray_from_data(self.npar, 1, self._base.drop_square)
384
def __set__(self, drop_sq):
387
d_ndr = numpy.atleast_1d(narray(drop_sq, copy=False, subok=True,
389
for i from 0 <= i < min(self.npar, d_ndr.size):
390
self._base.drop_square[i] = d_ndr[i]
392
def update(self, **modelargs):
393
family = modelargs.get('family', None)
394
if family is not None:
397
span = modelargs.get('span', None)
401
degree = modelargs.get('degree', None)
402
if degree is not None:
405
normalize = modelargs.get('normalize', None)
406
if normalize is not None:
407
self.normalize = normalize
409
parametric = modelargs.get('parametric', None)
410
if parametric is not None:
411
self.parametric_flags = parametric
413
drop_square = modelargs.get('drop_square', None)
414
if drop_square is not None:
415
self.drop_square_flags = drop_square
418
return "<loess object: model parameters>"
421
strg = ["Model parameters.....",
422
"Family : %s" % self.family,
423
"Span : %s" % self.span,
424
"Degree : %s" % self.degree,
425
"Normalized : %s" % self.normalize,
426
"Parametric : %s" % self.parametric_flags[:self.npar],
427
"Drop_square : %s" % self.drop_square_flags[:self.npar]
429
return '\n'.join(strg)
431
#####---------------------------------------------------------------------------
432
#---- ---- loess outputs ---
433
#####---------------------------------------------------------------------------
434
cdef class loess_outputs:
435
"""Outputs of a loess fit. This object is automatically created with empty
436
values when a new loess object is instantiated. The object gets filled when the
437
loess.fit() method is called.
440
fitted_values : ndarray
441
The (n,) ndarray of fitted values.
442
fitted_residuals : ndarray
443
The (n,) ndarray of fitted residuals (observations - fitted values).
445
Equivalent number of parameters.
447
Estimate of the scale of residuals.
449
Statistical parameter used in the computation of standard errors.
451
Statistical parameter used in the computation of standard errors.
452
pseudovalues : ndarray
453
The (n,) ndarray of adjusted values of the response when robust estimation
456
Trace of the operator hat matrix.
458
Diagonal of the operator hat matrix.
460
The (n,) ndarray of robustness weights for robust fitting.
462
The (p,) array of normalization divisors for numeric predictors.
464
cdef c_loess.c_loess_outputs *_base
466
cdef readonly int activated
467
cdef setup(self,c_loess.c_loess_outputs *base, long nobs, long npar):
471
self.activated = False
473
property fitted_values:
475
return floatarray_from_data(self.nobs, 1, self._base.fitted_values)
477
property fitted_residuals:
479
return floatarray_from_data(self.nobs, 1, self._base.fitted_residuals)
481
property pseudovalues:
483
return floatarray_from_data(self.nobs, 1, self._base.pseudovalues)
487
return floatarray_from_data(self.nobs, 1, self._base.diagonal)
491
return floatarray_from_data(self.nobs, 1, self._base.robust)
495
return floatarray_from_data(self.npar, 1, self._base.divisor)
499
return self._base.enp
507
return self._base.one_delta
511
return self._base.two_delta
515
return self._base.trace_hat
518
strg = ["Outputs................",
519
"Fitted values : %s\n" % self.fitted_values,
520
"Fitted residuals : %s\n" % self.fitted_residuals,
521
"Eqv. nb of parameters : %s" % self.enp,
522
"Residual error : %s" % self.s,
523
"Deltas : %s - %s" % (self.one_delta, self.two_delta),
524
"Normalization factors : %s" % self.divisor,]
525
return '\n'.join(strg)
529
#####---------------------------------------------------------------------------
530
#---- ---- loess confidence ---
531
#####---------------------------------------------------------------------------
532
cdef class conf_intervals:
533
"""Pointwise confidence intervals of a loess-predicted object:
539
Lower bounds of the confidence intervals.
541
Upper bounds of the confidence intervals.
543
cdef c_loess.c_conf_inv _base
544
cdef readonly ndarray lower, fit, upper
546
def __dealloc__(self):
547
c_loess.pw_free_mem(&self._base)
549
cdef setup(self, c_loess.c_conf_inv base, long nest):
551
self.fit = floatarray_from_data(nest, 1, base.fit)
552
self.upper = floatarray_from_data(nest, 1, base.upper)
553
self.lower = floatarray_from_data(nest, 1, base.lower)
557
tmp_ndr = numpy.r_[[self.lower,self.fit,self.upper]].T
558
return "Confidence intervals....\nLower b./ fit / upper b.\n%s" % \
561
#####---------------------------------------------------------------------------
562
#---- ---- loess predictions ---
563
#####---------------------------------------------------------------------------
564
cdef class loess_predicted:
565
"""Predicted values and standard errors of a loess object
569
The (m,) ndarray of loess values evaluated at newdata
571
The (m,) ndarray of the estimates of the standard error on the estimated
573
residual_scale : float
574
Estimate of the scale of the residuals
576
Degrees of freedom of the t-distribution used to compute pointwise
577
confidence intervals for the evaluated surface.
579
Number of new observations.
582
cdef c_loess.c_prediction _base
583
cdef readonly long nest
584
cdef readonly conf_intervals confidence_intervals
586
def __dealloc__(self):
587
c_loess.pred_free_mem(&self._base)
589
cdef setup(self, c_loess.c_prediction base, long nest):
595
return floatarray_from_data(self.nest, 1, self._base.fit)
599
return floatarray_from_data(self.nest, 1, self._base.se_fit)
601
property residual_scale:
603
return self._base.residual_scale
609
def confidence(self, coverage=0.95):
610
"""Returns the pointwise confidence intervals for each predicted values,
611
at the given confidence interval coverage.
615
Confidence level of the confidence intervals limits, as a fraction.
618
A new conf_intervals object, consisting of:
622
Lower bounds of the confidence intervals.
624
Upper bounds of the confidence intervals.
626
cdef c_loess.c_conf_inv _confintv
628
coverage = 1. - coverage
630
raise ValueError("The coverage precentage should be between 0 and 1!")
631
c_loess.c_pointwise(&self._base, self.nest, coverage, &_confintv)
632
self.confidence_intervals = conf_intervals()
633
self.confidence_intervals.setup(_confintv, self.nest)
634
return self.confidence_intervals
637
strg = ["Outputs................",
638
"Predicted values : %s\n" % self.values,
639
"Predicted std error : %s\n" % self.stderr,
640
"Residual scale : %s" % self.residual_scale,
641
"Degrees of freedom : %s" % self.df,
642
# "Confidence intervals : %s" % self.confidence,
644
return '\n'.join(strg)
647
#####---------------------------------------------------------------------------
648
#---- ---- loess base class ---
649
#####---------------------------------------------------------------------------
655
A (n,p) ndarray of independent variables, with n the number of observations
656
and p the number of variables.
658
A (n,) ndarray of observations
660
A (n,) ndarray of weights to be given to individual observations in the
661
sum of squared residuals that forms the local fitting criterion. If not
662
None, the weights should be non negative. If the different observations
663
have non-equal variances, the weights should be inversely proportional
665
By default, an unweighted fit is carried out (all the weights are one).
666
surface : string ["interpolate"]
667
Determines whether the fitted surface is computed directly at all points
668
("direct") or whether an interpolation method is used ("interpolate").
669
The default ("interpolate") is what most users should use unless special
670
circumstances warrant.
671
statistics : string ["approximate"]
672
Determines whether the statistical quantities are computed exactly
673
("exact") or approximately ("approximate"). "exact" should only be used
674
for testing the approximation in statistical development and is not meant
675
for routine usage because computation time can be horrendous.
676
trace_hat : string ["wait.to.decide"]
677
Determines how the trace of the hat matrix should be computed. The hat
678
matrix is used in the computation of the statistical quantities.
679
If "exact", an exact computation is done; this could be slow when the
680
number of observations n becomes large. If "wait.to.decide" is selected,
681
then a default is "exact" for n < 500 and "approximate" otherwise.
682
This option is only useful when the fitted surface is interpolated. If
683
surface is "exact", an exact computation is always done for the trace.
684
Setting trace_hat to "approximate" for large dataset will substantially
685
reduce the computation time.
687
Number of iterations of the robust fitting method. If the family is
688
"gaussian", the number of iterations is set to 0.
690
Maximum cell size of the kd-tree. Suppose k = floor(n*cell*span),
691
where n is the number of observations, and span the smoothing parameter.
692
Then, a cell is further divided if the number of observations within it
693
is greater than or equal to k. This option is only used if the surface
696
Smoothing factor, as a fraction of the number of points to take into
699
Overall degree of locally-fitted polynomial. 1 is locally-linear
700
fitting and 2 is locally-quadratic fitting. Degree should be 2 at most.
701
normalize : boolean [True]
702
Determines whether the independent variables should be normalized.
703
If True, the normalization is performed by setting the 10% trimmed
704
standard deviation to one. If False, no normalization is carried out.
705
This option is only useful for more than one variable. For spatial
706
coordinates predictors or variables with a common scale, it should be
708
family : string ["gaussian"]
709
Determines the assumed distribution of the errors. The values are
710
"gaussian" or "symmetric". If "gaussian" is selected, the fit is
711
performed with least-squares. If "symmetric" is selected, the fit
712
is performed robustly by redescending M-estimators.
713
parametric_flags : sequence [ [False]*p ]
714
Indicates which independent variables should be conditionally-parametric
715
(if there are two or more independent variables). The argument should
716
be a sequence of booleans, with the same size as the number of independent
717
variables, specified in the order of the predictor group ordered in x.
718
drop_square : sequence [ [False]* p]
719
When there are two or more independent variables and when a 2nd order
720
polynomial is used, "drop_square_flags" specifies those numeric predictors
721
whose squares should be dropped from the set of fitting variables.
722
The method of specification is the same as for parametric.
725
fitted_values : ndarray
726
The (n,) ndarray of fitted values.
727
fitted_residuals : ndarray
728
The (n,) ndarray of fitted residuals (observations - fitted values).
730
Equivalent number of parameters.
732
Estimate of the scale of residuals.
734
Statistical parameter used in the computation of standard errors.
736
Statistical parameter used in the computation of standard errors.
737
pseudovalues : ndarray
738
The (n,) ndarray of adjusted values of the response when robust estimation
741
Trace of the operator hat matrix.
743
Diagonal of the operator hat matrix.
745
The (n,) ndarray of robustness weights for robust fitting.
747
The (p,) array of normalization divisors for numeric predictors.
750
cdef c_loess.c_loess _base
751
cdef readonly loess_inputs inputs
752
cdef readonly loess_model model
753
cdef readonly loess_control control
754
cdef readonly loess_kd_tree kd_tree
755
cdef readonly loess_outputs outputs
756
cdef readonly loess_predicted predicted
758
def __init__(self, object x, object y, object weights=None, **options):
760
cdef ndarray x_ndr, y_ndr
761
cdef double *x_dat, *y_dat
763
# Initialize the inputs .......
764
self.inputs = loess_inputs(x, y)
765
x_dat = <double *>self.inputs.x_eff.data
766
y_dat = <double *>self.inputs.y_eff.data
769
c_loess.loess_setup(x_dat, y_dat, n, p, &self._base)
770
# Sets the _base input .........
771
self.inputs._base = &self._base.inputs
772
# Initialize the model parameters
773
self.model = loess_model()
774
self.model.setup(&self._base.model, p)
775
# Initialize the control parameters
776
self.control = loess_control()
777
self.control._base = &self._base.control
778
# Initialize the kd tree ......
779
self.kd_tree = loess_kd_tree()
780
self.kd_tree._base = &self._base.kd_tree
781
# Initialize the outputs ......
782
self.outputs = loess_outputs()
783
self.outputs.setup(&self._base.outputs, n, p,)
784
# Process options .............
787
for (k,v) in options.iteritems():
788
if k in ('family', 'span', 'degree', 'normalize',
789
'parametric', 'drop_square',):
791
elif k in ('surface', 'statistics', 'trace_hat',
792
'iterations', 'cell'):
794
self.control.update(**controlopt)
795
self.model.update(**modelopt)
796
#......................................................
798
"""Computes the loess parameters on the current inputs and sets of parameters."""
799
c_loess.loess_fit(&self._base)
800
self.outputs.activated = True
801
if self._base.status.err_status:
802
raise ValueError(self._base.status.err_msg)
804
#......................................................
805
def input_summary(self):
806
"""Returns some generic information about the loess parameters.
808
toprint = [str(self.inputs), str(self.model), str(self.control)]
809
return "\n".join(toprint)
811
def output_summary(self):
812
"""Returns some generic information about the loess fit."""
813
print "Number of Observations : %d" % self.inputs.nobs
814
print "Fit flag : %d" % bool(self.outputs.activated)
815
print "Equivalent Number of Parameters: %.1f" % self.outputs.enp
816
if self.model.family == "gaussian":
817
print "Residual Standard Error : %.4f" % self.outputs.s
819
print "Residual Scale Estimate : %.4f" % self.outputs.s
820
#......................................................
821
def predict(self, newdata, stderror=False):
822
"""Computes loess estimates at the given new data points newdata. Returns
823
a loess_predicted object, whose attributes are described below.
827
The (m,p) array of independent variables where the surface must be estimated,
828
with m the number of new data points, and p the number of independent
831
Whether the standard error should be computed
834
A new loess_predicted object, consisting of:
836
The (m,) ndarray of loess values evaluated at newdata
838
The (m,) ndarray of the estimates of the standard error on the estimated
840
residual_scale : float
841
Estimate of the scale of the residuals
843
Degrees of freedom of the t-distribution used to compute pointwise
844
confidence intervals for the evaluated surface.
846
Number of new observations.
850
cdef c_loess.c_prediction _prediction
852
# Make sure there's been a fit earlier ...
853
if self.outputs.activated == 0:
854
c_loess.loess_fit(&self._base)
855
self.outputs.activated = True
856
if self._base.status.err_status:
857
raise ValueError(self._base.status.err_msg)
858
# Note : we need a copy as we may have to normalize
859
p_ndr = narray(newdata, copy=True, subok=True, order='C').ravel()
860
p_dat = <double *>p_ndr.data
861
# Test the compatibility of sizes .......
863
raise ValueError("Can't predict without input data !")
864
(m, notOK) = divmod(len(p_ndr), self.inputs.npar)
867
"Incompatible data size: there should be as many rows as parameters")
869
c_loess.c_predict(p_dat, m, &self._base, &_prediction, stderror)
870
if self._base.status.err_status:
871
raise ValueError(self._base.status.err_msg)
872
self.predicted = loess_predicted()
873
self.predicted._base = _prediction
874
self.predicted.nest = m
875
# self.predicted.setup(_prediction, m)
876
return self.predicted
878
def __dealloc__(self):
879
c_loess.loess_free_mem(&self._base)
880
#......................................................
883
#####---------------------------------------------------------------------------
884
#---- ---- loess anova ---
885
#####---------------------------------------------------------------------------
887
cdef readonly double dfn, dfd, F_value, Pr_F
889
def __init__(self, loess_one, loess_two):
890
cdef double one_d1, one_d2, one_s, two_d1, two_d2, two_s, rssdiff,\
891
d1diff, tmp, df1, df2
893
if not isinstance(loess_one, loess) or not isinstance(loess_two, loess):
894
raise ValueError("Arguments should be valid loess objects!"\
895
"got '%s' instead" % type(loess_one))
897
out_one = loess_one.outputs
898
out_two = loess_two.outputs
900
one_d1 = out_one.one_delta
901
one_d2 = out_one.two_delta
904
two_d1 = out_two.one_delta
905
two_d2 = out_two.two_delta
908
rssdiff = abs(one_s * one_s * one_d1 - two_s * two_s * two_d1)
909
d1diff = abs(one_d1 - two_d1)
910
self.dfn = d1diff * d1diff / abs(one_d2 - two_d2)
913
if out_one.enp > out_two.enp:
914
self.dfd = one_d1 * one_d1 / one_d2
917
self.dfd = two_d1 * two_d1 / two_d2
920
F_value = (rssdiff / d1diff) / (tmp * tmp)
922
self.Pr_F = 1. - c_loess.ibeta(F_value*df1/(df2+F_value*df1), df1/2, df2/2)
923
self.F_value = F_value
b'\\ No newline at end of file'
added
removed
scipy/sandbox/pyloess/src/_loess.pyx
