comparison DEPENDENCIES/mingw32/Python27/Lib/site-packages/numpy/polynomial/hermite.py @ 87:2a2c65a20a8b

Add Python libs and headers
author Chris Cannam
date Wed, 25 Feb 2015 14:05:22 +0000
parents
children
comparison
equal deleted inserted replaced
86:413a9d26189e 87:2a2c65a20a8b
1 """
2 Objects for dealing with Hermite series.
3
4 This module provides a number of objects (mostly functions) useful for
5 dealing with Hermite series, including a `Hermite` class that
6 encapsulates the usual arithmetic operations. (General information
7 on how this module represents and works with such polynomials is in the
8 docstring for its "parent" sub-package, `numpy.polynomial`).
9
10 Constants
11 ---------
12 - `hermdomain` -- Hermite series default domain, [-1,1].
13 - `hermzero` -- Hermite series that evaluates identically to 0.
14 - `hermone` -- Hermite series that evaluates identically to 1.
15 - `hermx` -- Hermite series for the identity map, ``f(x) = x``.
16
17 Arithmetic
18 ----------
19 - `hermmulx` -- multiply a Hermite series in ``P_i(x)`` by ``x``.
20 - `hermadd` -- add two Hermite series.
21 - `hermsub` -- subtract one Hermite series from another.
22 - `hermmul` -- multiply two Hermite series.
23 - `hermdiv` -- divide one Hermite series by another.
24 - `hermval` -- evaluate a Hermite series at given points.
25 - `hermval2d` -- evaluate a 2D Hermite series at given points.
26 - `hermval3d` -- evaluate a 3D Hermite series at given points.
27 - `hermgrid2d` -- evaluate a 2D Hermite series on a Cartesian product.
28 - `hermgrid3d` -- evaluate a 3D Hermite series on a Cartesian product.
29
30 Calculus
31 --------
32 - `hermder` -- differentiate a Hermite series.
33 - `hermint` -- integrate a Hermite series.
34
35 Misc Functions
36 --------------
37 - `hermfromroots` -- create a Hermite series with specified roots.
38 - `hermroots` -- find the roots of a Hermite series.
39 - `hermvander` -- Vandermonde-like matrix for Hermite polynomials.
40 - `hermvander2d` -- Vandermonde-like matrix for 2D power series.
41 - `hermvander3d` -- Vandermonde-like matrix for 3D power series.
42 - `hermgauss` -- Gauss-Hermite quadrature, points and weights.
43 - `hermweight` -- Hermite weight function.
44 - `hermcompanion` -- symmetrized companion matrix in Hermite form.
45 - `hermfit` -- least-squares fit returning a Hermite series.
46 - `hermtrim` -- trim leading coefficients from a Hermite series.
47 - `hermline` -- Hermite series of given straight line.
48 - `herm2poly` -- convert a Hermite series to a polynomial.
49 - `poly2herm` -- convert a polynomial to a Hermite series.
50
51 Classes
52 -------
53 - `Hermite` -- A Hermite series class.
54
55 See also
56 --------
57 `numpy.polynomial`
58
59 """
60 from __future__ import division, absolute_import, print_function
61
62 import warnings
63 import numpy as np
64 import numpy.linalg as la
65
66 from . import polyutils as pu
67 from ._polybase import ABCPolyBase
68
69 __all__ = [
70 'hermzero', 'hermone', 'hermx', 'hermdomain', 'hermline', 'hermadd',
71 'hermsub', 'hermmulx', 'hermmul', 'hermdiv', 'hermpow', 'hermval',
72 'hermder', 'hermint', 'herm2poly', 'poly2herm', 'hermfromroots',
73 'hermvander', 'hermfit', 'hermtrim', 'hermroots', 'Hermite',
74 'hermval2d', 'hermval3d', 'hermgrid2d', 'hermgrid3d', 'hermvander2d',
75 'hermvander3d', 'hermcompanion', 'hermgauss', 'hermweight']
76
77 hermtrim = pu.trimcoef
78
79
80 def poly2herm(pol):
81 """
82 poly2herm(pol)
83
84 Convert a polynomial to a Hermite series.
85
86 Convert an array representing the coefficients of a polynomial (relative
87 to the "standard" basis) ordered from lowest degree to highest, to an
88 array of the coefficients of the equivalent Hermite series, ordered
89 from lowest to highest degree.
90
91 Parameters
92 ----------
93 pol : array_like
94 1-D array containing the polynomial coefficients
95
96 Returns
97 -------
98 c : ndarray
99 1-D array containing the coefficients of the equivalent Hermite
100 series.
101
102 See Also
103 --------
104 herm2poly
105
106 Notes
107 -----
108 The easy way to do conversions between polynomial basis sets
109 is to use the convert method of a class instance.
110
111 Examples
112 --------
113 >>> from numpy.polynomial.hermite import poly2herm
114 >>> poly2herm(np.arange(4))
115 array([ 1. , 2.75 , 0.5 , 0.375])
116
117 """
118 [pol] = pu.as_series([pol])
119 deg = len(pol) - 1
120 res = 0
121 for i in range(deg, -1, -1):
122 res = hermadd(hermmulx(res), pol[i])
123 return res
124
125
126 def herm2poly(c):
127 """
128 Convert a Hermite series to a polynomial.
129
130 Convert an array representing the coefficients of a Hermite series,
131 ordered from lowest degree to highest, to an array of the coefficients
132 of the equivalent polynomial (relative to the "standard" basis) ordered
133 from lowest to highest degree.
134
135 Parameters
136 ----------
137 c : array_like
138 1-D array containing the Hermite series coefficients, ordered
139 from lowest order term to highest.
140
141 Returns
142 -------
143 pol : ndarray
144 1-D array containing the coefficients of the equivalent polynomial
145 (relative to the "standard" basis) ordered from lowest order term
146 to highest.
147
148 See Also
149 --------
150 poly2herm
151
152 Notes
153 -----
154 The easy way to do conversions between polynomial basis sets
155 is to use the convert method of a class instance.
156
157 Examples
158 --------
159 >>> from numpy.polynomial.hermite import herm2poly
160 >>> herm2poly([ 1. , 2.75 , 0.5 , 0.375])
161 array([ 0., 1., 2., 3.])
162
163 """
164 from .polynomial import polyadd, polysub, polymulx
165
166 [c] = pu.as_series([c])
167 n = len(c)
168 if n == 1:
169 return c
170 if n == 2:
171 c[1] *= 2
172 return c
173 else:
174 c0 = c[-2]
175 c1 = c[-1]
176 # i is the current degree of c1
177 for i in range(n - 1, 1, -1):
178 tmp = c0
179 c0 = polysub(c[i - 2], c1*(2*(i - 1)))
180 c1 = polyadd(tmp, polymulx(c1)*2)
181 return polyadd(c0, polymulx(c1)*2)
182
183 #
184 # These are constant arrays are of integer type so as to be compatible
185 # with the widest range of other types, such as Decimal.
186 #
187
188 # Hermite
189 hermdomain = np.array([-1, 1])
190
191 # Hermite coefficients representing zero.
192 hermzero = np.array([0])
193
194 # Hermite coefficients representing one.
195 hermone = np.array([1])
196
197 # Hermite coefficients representing the identity x.
198 hermx = np.array([0, 1/2])
199
200
201 def hermline(off, scl):
202 """
203 Hermite series whose graph is a straight line.
204
205
206
207 Parameters
208 ----------
209 off, scl : scalars
210 The specified line is given by ``off + scl*x``.
211
212 Returns
213 -------
214 y : ndarray
215 This module's representation of the Hermite series for
216 ``off + scl*x``.
217
218 See Also
219 --------
220 polyline, chebline
221
222 Examples
223 --------
224 >>> from numpy.polynomial.hermite import hermline, hermval
225 >>> hermval(0,hermline(3, 2))
226 3.0
227 >>> hermval(1,hermline(3, 2))
228 5.0
229
230 """
231 if scl != 0:
232 return np.array([off, scl/2])
233 else:
234 return np.array([off])
235
236
237 def hermfromroots(roots):
238 """
239 Generate a Hermite series with given roots.
240
241 The function returns the coefficients of the polynomial
242
243 .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
244
245 in Hermite form, where the `r_n` are the roots specified in `roots`.
246 If a zero has multiplicity n, then it must appear in `roots` n times.
247 For instance, if 2 is a root of multiplicity three and 3 is a root of
248 multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The
249 roots can appear in any order.
250
251 If the returned coefficients are `c`, then
252
253 .. math:: p(x) = c_0 + c_1 * H_1(x) + ... + c_n * H_n(x)
254
255 The coefficient of the last term is not generally 1 for monic
256 polynomials in Hermite form.
257
258 Parameters
259 ----------
260 roots : array_like
261 Sequence containing the roots.
262
263 Returns
264 -------
265 out : ndarray
266 1-D array of coefficients. If all roots are real then `out` is a
267 real array, if some of the roots are complex, then `out` is complex
268 even if all the coefficients in the result are real (see Examples
269 below).
270
271 See Also
272 --------
273 polyfromroots, legfromroots, lagfromroots, chebfromroots,
274 hermefromroots.
275
276 Examples
277 --------
278 >>> from numpy.polynomial.hermite import hermfromroots, hermval
279 >>> coef = hermfromroots((-1, 0, 1))
280 >>> hermval((-1, 0, 1), coef)
281 array([ 0., 0., 0.])
282 >>> coef = hermfromroots((-1j, 1j))
283 >>> hermval((-1j, 1j), coef)
284 array([ 0.+0.j, 0.+0.j])
285
286 """
287 if len(roots) == 0:
288 return np.ones(1)
289 else:
290 [roots] = pu.as_series([roots], trim=False)
291 roots.sort()
292 p = [hermline(-r, 1) for r in roots]
293 n = len(p)
294 while n > 1:
295 m, r = divmod(n, 2)
296 tmp = [hermmul(p[i], p[i+m]) for i in range(m)]
297 if r:
298 tmp[0] = hermmul(tmp[0], p[-1])
299 p = tmp
300 n = m
301 return p[0]
302
303
304 def hermadd(c1, c2):
305 """
306 Add one Hermite series to another.
307
308 Returns the sum of two Hermite series `c1` + `c2`. The arguments
309 are sequences of coefficients ordered from lowest order term to
310 highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
311
312 Parameters
313 ----------
314 c1, c2 : array_like
315 1-D arrays of Hermite series coefficients ordered from low to
316 high.
317
318 Returns
319 -------
320 out : ndarray
321 Array representing the Hermite series of their sum.
322
323 See Also
324 --------
325 hermsub, hermmul, hermdiv, hermpow
326
327 Notes
328 -----
329 Unlike multiplication, division, etc., the sum of two Hermite series
330 is a Hermite series (without having to "reproject" the result onto
331 the basis set) so addition, just like that of "standard" polynomials,
332 is simply "component-wise."
333
334 Examples
335 --------
336 >>> from numpy.polynomial.hermite import hermadd
337 >>> hermadd([1, 2, 3], [1, 2, 3, 4])
338 array([ 2., 4., 6., 4.])
339
340 """
341 # c1, c2 are trimmed copies
342 [c1, c2] = pu.as_series([c1, c2])
343 if len(c1) > len(c2):
344 c1[:c2.size] += c2
345 ret = c1
346 else:
347 c2[:c1.size] += c1
348 ret = c2
349 return pu.trimseq(ret)
350
351
352 def hermsub(c1, c2):
353 """
354 Subtract one Hermite series from another.
355
356 Returns the difference of two Hermite series `c1` - `c2`. The
357 sequences of coefficients are from lowest order term to highest, i.e.,
358 [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
359
360 Parameters
361 ----------
362 c1, c2 : array_like
363 1-D arrays of Hermite series coefficients ordered from low to
364 high.
365
366 Returns
367 -------
368 out : ndarray
369 Of Hermite series coefficients representing their difference.
370
371 See Also
372 --------
373 hermadd, hermmul, hermdiv, hermpow
374
375 Notes
376 -----
377 Unlike multiplication, division, etc., the difference of two Hermite
378 series is a Hermite series (without having to "reproject" the result
379 onto the basis set) so subtraction, just like that of "standard"
380 polynomials, is simply "component-wise."
381
382 Examples
383 --------
384 >>> from numpy.polynomial.hermite import hermsub
385 >>> hermsub([1, 2, 3, 4], [1, 2, 3])
386 array([ 0., 0., 0., 4.])
387
388 """
389 # c1, c2 are trimmed copies
390 [c1, c2] = pu.as_series([c1, c2])
391 if len(c1) > len(c2):
392 c1[:c2.size] -= c2
393 ret = c1
394 else:
395 c2 = -c2
396 c2[:c1.size] += c1
397 ret = c2
398 return pu.trimseq(ret)
399
400
401 def hermmulx(c):
402 """Multiply a Hermite series by x.
403
404 Multiply the Hermite series `c` by x, where x is the independent
405 variable.
406
407
408 Parameters
409 ----------
410 c : array_like
411 1-D array of Hermite series coefficients ordered from low to
412 high.
413
414 Returns
415 -------
416 out : ndarray
417 Array representing the result of the multiplication.
418
419 Notes
420 -----
421 The multiplication uses the recursion relationship for Hermite
422 polynomials in the form
423
424 .. math::
425
426 xP_i(x) = (P_{i + 1}(x)/2 + i*P_{i - 1}(x))
427
428 Examples
429 --------
430 >>> from numpy.polynomial.hermite import hermmulx
431 >>> hermmulx([1, 2, 3])
432 array([ 2. , 6.5, 1. , 1.5])
433
434 """
435 # c is a trimmed copy
436 [c] = pu.as_series([c])
437 # The zero series needs special treatment
438 if len(c) == 1 and c[0] == 0:
439 return c
440
441 prd = np.empty(len(c) + 1, dtype=c.dtype)
442 prd[0] = c[0]*0
443 prd[1] = c[0]/2
444 for i in range(1, len(c)):
445 prd[i + 1] = c[i]/2
446 prd[i - 1] += c[i]*i
447 return prd
448
449
450 def hermmul(c1, c2):
451 """
452 Multiply one Hermite series by another.
453
454 Returns the product of two Hermite series `c1` * `c2`. The arguments
455 are sequences of coefficients, from lowest order "term" to highest,
456 e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
457
458 Parameters
459 ----------
460 c1, c2 : array_like
461 1-D arrays of Hermite series coefficients ordered from low to
462 high.
463
464 Returns
465 -------
466 out : ndarray
467 Of Hermite series coefficients representing their product.
468
469 See Also
470 --------
471 hermadd, hermsub, hermdiv, hermpow
472
473 Notes
474 -----
475 In general, the (polynomial) product of two C-series results in terms
476 that are not in the Hermite polynomial basis set. Thus, to express
477 the product as a Hermite series, it is necessary to "reproject" the
478 product onto said basis set, which may produce "unintuitive" (but
479 correct) results; see Examples section below.
480
481 Examples
482 --------
483 >>> from numpy.polynomial.hermite import hermmul
484 >>> hermmul([1, 2, 3], [0, 1, 2])
485 array([ 52., 29., 52., 7., 6.])
486
487 """
488 # s1, s2 are trimmed copies
489 [c1, c2] = pu.as_series([c1, c2])
490
491 if len(c1) > len(c2):
492 c = c2
493 xs = c1
494 else:
495 c = c1
496 xs = c2
497
498 if len(c) == 1:
499 c0 = c[0]*xs
500 c1 = 0
501 elif len(c) == 2:
502 c0 = c[0]*xs
503 c1 = c[1]*xs
504 else:
505 nd = len(c)
506 c0 = c[-2]*xs
507 c1 = c[-1]*xs
508 for i in range(3, len(c) + 1):
509 tmp = c0
510 nd = nd - 1
511 c0 = hermsub(c[-i]*xs, c1*(2*(nd - 1)))
512 c1 = hermadd(tmp, hermmulx(c1)*2)
513 return hermadd(c0, hermmulx(c1)*2)
514
515
516 def hermdiv(c1, c2):
517 """
518 Divide one Hermite series by another.
519
520 Returns the quotient-with-remainder of two Hermite series
521 `c1` / `c2`. The arguments are sequences of coefficients from lowest
522 order "term" to highest, e.g., [1,2,3] represents the series
523 ``P_0 + 2*P_1 + 3*P_2``.
524
525 Parameters
526 ----------
527 c1, c2 : array_like
528 1-D arrays of Hermite series coefficients ordered from low to
529 high.
530
531 Returns
532 -------
533 [quo, rem] : ndarrays
534 Of Hermite series coefficients representing the quotient and
535 remainder.
536
537 See Also
538 --------
539 hermadd, hermsub, hermmul, hermpow
540
541 Notes
542 -----
543 In general, the (polynomial) division of one Hermite series by another
544 results in quotient and remainder terms that are not in the Hermite
545 polynomial basis set. Thus, to express these results as a Hermite
546 series, it is necessary to "reproject" the results onto the Hermite
547 basis set, which may produce "unintuitive" (but correct) results; see
548 Examples section below.
549
550 Examples
551 --------
552 >>> from numpy.polynomial.hermite import hermdiv
553 >>> hermdiv([ 52., 29., 52., 7., 6.], [0, 1, 2])
554 (array([ 1., 2., 3.]), array([ 0.]))
555 >>> hermdiv([ 54., 31., 52., 7., 6.], [0, 1, 2])
556 (array([ 1., 2., 3.]), array([ 2., 2.]))
557 >>> hermdiv([ 53., 30., 52., 7., 6.], [0, 1, 2])
558 (array([ 1., 2., 3.]), array([ 1., 1.]))
559
560 """
561 # c1, c2 are trimmed copies
562 [c1, c2] = pu.as_series([c1, c2])
563 if c2[-1] == 0:
564 raise ZeroDivisionError()
565
566 lc1 = len(c1)
567 lc2 = len(c2)
568 if lc1 < lc2:
569 return c1[:1]*0, c1
570 elif lc2 == 1:
571 return c1/c2[-1], c1[:1]*0
572 else:
573 quo = np.empty(lc1 - lc2 + 1, dtype=c1.dtype)
574 rem = c1
575 for i in range(lc1 - lc2, - 1, -1):
576 p = hermmul([0]*i + [1], c2)
577 q = rem[-1]/p[-1]
578 rem = rem[:-1] - q*p[:-1]
579 quo[i] = q
580 return quo, pu.trimseq(rem)
581
582
583 def hermpow(c, pow, maxpower=16):
584 """Raise a Hermite series to a power.
585
586 Returns the Hermite series `c` raised to the power `pow`. The
587 argument `c` is a sequence of coefficients ordered from low to high.
588 i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.``
589
590 Parameters
591 ----------
592 c : array_like
593 1-D array of Hermite series coefficients ordered from low to
594 high.
595 pow : integer
596 Power to which the series will be raised
597 maxpower : integer, optional
598 Maximum power allowed. This is mainly to limit growth of the series
599 to unmanageable size. Default is 16
600
601 Returns
602 -------
603 coef : ndarray
604 Hermite series of power.
605
606 See Also
607 --------
608 hermadd, hermsub, hermmul, hermdiv
609
610 Examples
611 --------
612 >>> from numpy.polynomial.hermite import hermpow
613 >>> hermpow([1, 2, 3], 2)
614 array([ 81., 52., 82., 12., 9.])
615
616 """
617 # c is a trimmed copy
618 [c] = pu.as_series([c])
619 power = int(pow)
620 if power != pow or power < 0:
621 raise ValueError("Power must be a non-negative integer.")
622 elif maxpower is not None and power > maxpower:
623 raise ValueError("Power is too large")
624 elif power == 0:
625 return np.array([1], dtype=c.dtype)
626 elif power == 1:
627 return c
628 else:
629 # This can be made more efficient by using powers of two
630 # in the usual way.
631 prd = c
632 for i in range(2, power + 1):
633 prd = hermmul(prd, c)
634 return prd
635
636
637 def hermder(c, m=1, scl=1, axis=0):
638 """
639 Differentiate a Hermite series.
640
641 Returns the Hermite series coefficients `c` differentiated `m` times
642 along `axis`. At each iteration the result is multiplied by `scl` (the
643 scaling factor is for use in a linear change of variable). The argument
644 `c` is an array of coefficients from low to high degree along each
645 axis, e.g., [1,2,3] represents the series ``1*H_0 + 2*H_1 + 3*H_2``
646 while [[1,2],[1,2]] represents ``1*H_0(x)*H_0(y) + 1*H_1(x)*H_0(y) +
647 2*H_0(x)*H_1(y) + 2*H_1(x)*H_1(y)`` if axis=0 is ``x`` and axis=1 is
648 ``y``.
649
650 Parameters
651 ----------
652 c : array_like
653 Array of Hermite series coefficients. If `c` is multidimensional the
654 different axis correspond to different variables with the degree in
655 each axis given by the corresponding index.
656 m : int, optional
657 Number of derivatives taken, must be non-negative. (Default: 1)
658 scl : scalar, optional
659 Each differentiation is multiplied by `scl`. The end result is
660 multiplication by ``scl**m``. This is for use in a linear change of
661 variable. (Default: 1)
662 axis : int, optional
663 Axis over which the derivative is taken. (Default: 0).
664
665 .. versionadded:: 1.7.0
666
667 Returns
668 -------
669 der : ndarray
670 Hermite series of the derivative.
671
672 See Also
673 --------
674 hermint
675
676 Notes
677 -----
678 In general, the result of differentiating a Hermite series does not
679 resemble the same operation on a power series. Thus the result of this
680 function may be "unintuitive," albeit correct; see Examples section
681 below.
682
683 Examples
684 --------
685 >>> from numpy.polynomial.hermite import hermder
686 >>> hermder([ 1. , 0.5, 0.5, 0.5])
687 array([ 1., 2., 3.])
688 >>> hermder([-0.5, 1./2., 1./8., 1./12., 1./16.], m=2)
689 array([ 1., 2., 3.])
690
691 """
692 c = np.array(c, ndmin=1, copy=1)
693 if c.dtype.char in '?bBhHiIlLqQpP':
694 c = c.astype(np.double)
695 cnt, iaxis = [int(t) for t in [m, axis]]
696
697 if cnt != m:
698 raise ValueError("The order of derivation must be integer")
699 if cnt < 0:
700 raise ValueError("The order of derivation must be non-negative")
701 if iaxis != axis:
702 raise ValueError("The axis must be integer")
703 if not -c.ndim <= iaxis < c.ndim:
704 raise ValueError("The axis is out of range")
705 if iaxis < 0:
706 iaxis += c.ndim
707
708 if cnt == 0:
709 return c
710
711 c = np.rollaxis(c, iaxis)
712 n = len(c)
713 if cnt >= n:
714 c = c[:1]*0
715 else:
716 for i in range(cnt):
717 n = n - 1
718 c *= scl
719 der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
720 for j in range(n, 0, -1):
721 der[j - 1] = (2*j)*c[j]
722 c = der
723 c = np.rollaxis(c, 0, iaxis + 1)
724 return c
725
726
727 def hermint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
728 """
729 Integrate a Hermite series.
730
731 Returns the Hermite series coefficients `c` integrated `m` times from
732 `lbnd` along `axis`. At each iteration the resulting series is
733 **multiplied** by `scl` and an integration constant, `k`, is added.
734 The scaling factor is for use in a linear change of variable. ("Buyer
735 beware": note that, depending on what one is doing, one may want `scl`
736 to be the reciprocal of what one might expect; for more information,
737 see the Notes section below.) The argument `c` is an array of
738 coefficients from low to high degree along each axis, e.g., [1,2,3]
739 represents the series ``H_0 + 2*H_1 + 3*H_2`` while [[1,2],[1,2]]
740 represents ``1*H_0(x)*H_0(y) + 1*H_1(x)*H_0(y) + 2*H_0(x)*H_1(y) +
741 2*H_1(x)*H_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
742
743 Parameters
744 ----------
745 c : array_like
746 Array of Hermite series coefficients. If c is multidimensional the
747 different axis correspond to different variables with the degree in
748 each axis given by the corresponding index.
749 m : int, optional
750 Order of integration, must be positive. (Default: 1)
751 k : {[], list, scalar}, optional
752 Integration constant(s). The value of the first integral at
753 ``lbnd`` is the first value in the list, the value of the second
754 integral at ``lbnd`` is the second value, etc. If ``k == []`` (the
755 default), all constants are set to zero. If ``m == 1``, a single
756 scalar can be given instead of a list.
757 lbnd : scalar, optional
758 The lower bound of the integral. (Default: 0)
759 scl : scalar, optional
760 Following each integration the result is *multiplied* by `scl`
761 before the integration constant is added. (Default: 1)
762 axis : int, optional
763 Axis over which the integral is taken. (Default: 0).
764
765 .. versionadded:: 1.7.0
766
767 Returns
768 -------
769 S : ndarray
770 Hermite series coefficients of the integral.
771
772 Raises
773 ------
774 ValueError
775 If ``m < 0``, ``len(k) > m``, ``np.isscalar(lbnd) == False``, or
776 ``np.isscalar(scl) == False``.
777
778 See Also
779 --------
780 hermder
781
782 Notes
783 -----
784 Note that the result of each integration is *multiplied* by `scl`.
785 Why is this important to note? Say one is making a linear change of
786 variable :math:`u = ax + b` in an integral relative to `x`. Then
787 .. math::`dx = du/a`, so one will need to set `scl` equal to
788 :math:`1/a` - perhaps not what one would have first thought.
789
790 Also note that, in general, the result of integrating a C-series needs
791 to be "reprojected" onto the C-series basis set. Thus, typically,
792 the result of this function is "unintuitive," albeit correct; see
793 Examples section below.
794
795 Examples
796 --------
797 >>> from numpy.polynomial.hermite import hermint
798 >>> hermint([1,2,3]) # integrate once, value 0 at 0.
799 array([ 1. , 0.5, 0.5, 0.5])
800 >>> hermint([1,2,3], m=2) # integrate twice, value & deriv 0 at 0
801 array([-0.5 , 0.5 , 0.125 , 0.08333333, 0.0625 ])
802 >>> hermint([1,2,3], k=1) # integrate once, value 1 at 0.
803 array([ 2. , 0.5, 0.5, 0.5])
804 >>> hermint([1,2,3], lbnd=-1) # integrate once, value 0 at -1
805 array([-2. , 0.5, 0.5, 0.5])
806 >>> hermint([1,2,3], m=2, k=[1,2], lbnd=-1)
807 array([ 1.66666667, -0.5 , 0.125 , 0.08333333, 0.0625 ])
808
809 """
810 c = np.array(c, ndmin=1, copy=1)
811 if c.dtype.char in '?bBhHiIlLqQpP':
812 c = c.astype(np.double)
813 if not np.iterable(k):
814 k = [k]
815 cnt, iaxis = [int(t) for t in [m, axis]]
816
817 if cnt != m:
818 raise ValueError("The order of integration must be integer")
819 if cnt < 0:
820 raise ValueError("The order of integration must be non-negative")
821 if len(k) > cnt:
822 raise ValueError("Too many integration constants")
823 if iaxis != axis:
824 raise ValueError("The axis must be integer")
825 if not -c.ndim <= iaxis < c.ndim:
826 raise ValueError("The axis is out of range")
827 if iaxis < 0:
828 iaxis += c.ndim
829
830 if cnt == 0:
831 return c
832
833 c = np.rollaxis(c, iaxis)
834 k = list(k) + [0]*(cnt - len(k))
835 for i in range(cnt):
836 n = len(c)
837 c *= scl
838 if n == 1 and np.all(c[0] == 0):
839 c[0] += k[i]
840 else:
841 tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
842 tmp[0] = c[0]*0
843 tmp[1] = c[0]/2
844 for j in range(1, n):
845 tmp[j + 1] = c[j]/(2*(j + 1))
846 tmp[0] += k[i] - hermval(lbnd, tmp)
847 c = tmp
848 c = np.rollaxis(c, 0, iaxis + 1)
849 return c
850
851
852 def hermval(x, c, tensor=True):
853 """
854 Evaluate an Hermite series at points x.
855
856 If `c` is of length `n + 1`, this function returns the value:
857
858 .. math:: p(x) = c_0 * H_0(x) + c_1 * H_1(x) + ... + c_n * H_n(x)
859
860 The parameter `x` is converted to an array only if it is a tuple or a
861 list, otherwise it is treated as a scalar. In either case, either `x`
862 or its elements must support multiplication and addition both with
863 themselves and with the elements of `c`.
864
865 If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If
866 `c` is multidimensional, then the shape of the result depends on the
867 value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
868 x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
869 scalars have shape (,).
870
871 Trailing zeros in the coefficients will be used in the evaluation, so
872 they should be avoided if efficiency is a concern.
873
874 Parameters
875 ----------
876 x : array_like, compatible object
877 If `x` is a list or tuple, it is converted to an ndarray, otherwise
878 it is left unchanged and treated as a scalar. In either case, `x`
879 or its elements must support addition and multiplication with
880 with themselves and with the elements of `c`.
881 c : array_like
882 Array of coefficients ordered so that the coefficients for terms of
883 degree n are contained in c[n]. If `c` is multidimensional the
884 remaining indices enumerate multiple polynomials. In the two
885 dimensional case the coefficients may be thought of as stored in
886 the columns of `c`.
887 tensor : boolean, optional
888 If True, the shape of the coefficient array is extended with ones
889 on the right, one for each dimension of `x`. Scalars have dimension 0
890 for this action. The result is that every column of coefficients in
891 `c` is evaluated for every element of `x`. If False, `x` is broadcast
892 over the columns of `c` for the evaluation. This keyword is useful
893 when `c` is multidimensional. The default value is True.
894
895 .. versionadded:: 1.7.0
896
897 Returns
898 -------
899 values : ndarray, algebra_like
900 The shape of the return value is described above.
901
902 See Also
903 --------
904 hermval2d, hermgrid2d, hermval3d, hermgrid3d
905
906 Notes
907 -----
908 The evaluation uses Clenshaw recursion, aka synthetic division.
909
910 Examples
911 --------
912 >>> from numpy.polynomial.hermite import hermval
913 >>> coef = [1,2,3]
914 >>> hermval(1, coef)
915 11.0
916 >>> hermval([[1,2],[3,4]], coef)
917 array([[ 11., 51.],
918 [ 115., 203.]])
919
920 """
921 c = np.array(c, ndmin=1, copy=0)
922 if c.dtype.char in '?bBhHiIlLqQpP':
923 c = c.astype(np.double)
924 if isinstance(x, (tuple, list)):
925 x = np.asarray(x)
926 if isinstance(x, np.ndarray) and tensor:
927 c = c.reshape(c.shape + (1,)*x.ndim)
928
929 x2 = x*2
930 if len(c) == 1:
931 c0 = c[0]
932 c1 = 0
933 elif len(c) == 2:
934 c0 = c[0]
935 c1 = c[1]
936 else:
937 nd = len(c)
938 c0 = c[-2]
939 c1 = c[-1]
940 for i in range(3, len(c) + 1):
941 tmp = c0
942 nd = nd - 1
943 c0 = c[-i] - c1*(2*(nd - 1))
944 c1 = tmp + c1*x2
945 return c0 + c1*x2
946
947
948 def hermval2d(x, y, c):
949 """
950 Evaluate a 2-D Hermite series at points (x, y).
951
952 This function returns the values:
953
954 .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * H_i(x) * H_j(y)
955
956 The parameters `x` and `y` are converted to arrays only if they are
957 tuples or a lists, otherwise they are treated as a scalars and they
958 must have the same shape after conversion. In either case, either `x`
959 and `y` or their elements must support multiplication and addition both
960 with themselves and with the elements of `c`.
961
962 If `c` is a 1-D array a one is implicitly appended to its shape to make
963 it 2-D. The shape of the result will be c.shape[2:] + x.shape.
964
965 Parameters
966 ----------
967 x, y : array_like, compatible objects
968 The two dimensional series is evaluated at the points `(x, y)`,
969 where `x` and `y` must have the same shape. If `x` or `y` is a list
970 or tuple, it is first converted to an ndarray, otherwise it is left
971 unchanged and if it isn't an ndarray it is treated as a scalar.
972 c : array_like
973 Array of coefficients ordered so that the coefficient of the term
974 of multi-degree i,j is contained in ``c[i,j]``. If `c` has
975 dimension greater than two the remaining indices enumerate multiple
976 sets of coefficients.
977
978 Returns
979 -------
980 values : ndarray, compatible object
981 The values of the two dimensional polynomial at points formed with
982 pairs of corresponding values from `x` and `y`.
983
984 See Also
985 --------
986 hermval, hermgrid2d, hermval3d, hermgrid3d
987
988 Notes
989 -----
990
991 .. versionadded::1.7.0
992
993 """
994 try:
995 x, y = np.array((x, y), copy=0)
996 except:
997 raise ValueError('x, y are incompatible')
998
999 c = hermval(x, c)
1000 c = hermval(y, c, tensor=False)
1001 return c
1002
1003
1004 def hermgrid2d(x, y, c):
1005 """
1006 Evaluate a 2-D Hermite series on the Cartesian product of x and y.
1007
1008 This function returns the values:
1009
1010 .. math:: p(a,b) = \sum_{i,j} c_{i,j} * H_i(a) * H_j(b)
1011
1012 where the points `(a, b)` consist of all pairs formed by taking
1013 `a` from `x` and `b` from `y`. The resulting points form a grid with
1014 `x` in the first dimension and `y` in the second.
1015
1016 The parameters `x` and `y` are converted to arrays only if they are
1017 tuples or a lists, otherwise they are treated as a scalars. In either
1018 case, either `x` and `y` or their elements must support multiplication
1019 and addition both with themselves and with the elements of `c`.
1020
1021 If `c` has fewer than two dimensions, ones are implicitly appended to
1022 its shape to make it 2-D. The shape of the result will be c.shape[2:] +
1023 x.shape.
1024
1025 Parameters
1026 ----------
1027 x, y : array_like, compatible objects
1028 The two dimensional series is evaluated at the points in the
1029 Cartesian product of `x` and `y`. If `x` or `y` is a list or
1030 tuple, it is first converted to an ndarray, otherwise it is left
1031 unchanged and, if it isn't an ndarray, it is treated as a scalar.
1032 c : array_like
1033 Array of coefficients ordered so that the coefficients for terms of
1034 degree i,j are contained in ``c[i,j]``. If `c` has dimension
1035 greater than two the remaining indices enumerate multiple sets of
1036 coefficients.
1037
1038 Returns
1039 -------
1040 values : ndarray, compatible object
1041 The values of the two dimensional polynomial at points in the Cartesian
1042 product of `x` and `y`.
1043
1044 See Also
1045 --------
1046 hermval, hermval2d, hermval3d, hermgrid3d
1047
1048 Notes
1049 -----
1050
1051 .. versionadded::1.7.0
1052
1053 """
1054 c = hermval(x, c)
1055 c = hermval(y, c)
1056 return c
1057
1058
1059 def hermval3d(x, y, z, c):
1060 """
1061 Evaluate a 3-D Hermite series at points (x, y, z).
1062
1063 This function returns the values:
1064
1065 .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * H_i(x) * H_j(y) * H_k(z)
1066
1067 The parameters `x`, `y`, and `z` are converted to arrays only if
1068 they are tuples or a lists, otherwise they are treated as a scalars and
1069 they must have the same shape after conversion. In either case, either
1070 `x`, `y`, and `z` or their elements must support multiplication and
1071 addition both with themselves and with the elements of `c`.
1072
1073 If `c` has fewer than 3 dimensions, ones are implicitly appended to its
1074 shape to make it 3-D. The shape of the result will be c.shape[3:] +
1075 x.shape.
1076
1077 Parameters
1078 ----------
1079 x, y, z : array_like, compatible object
1080 The three dimensional series is evaluated at the points
1081 `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If
1082 any of `x`, `y`, or `z` is a list or tuple, it is first converted
1083 to an ndarray, otherwise it is left unchanged and if it isn't an
1084 ndarray it is treated as a scalar.
1085 c : array_like
1086 Array of coefficients ordered so that the coefficient of the term of
1087 multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
1088 greater than 3 the remaining indices enumerate multiple sets of
1089 coefficients.
1090
1091 Returns
1092 -------
1093 values : ndarray, compatible object
1094 The values of the multidimensional polynomial on points formed with
1095 triples of corresponding values from `x`, `y`, and `z`.
1096
1097 See Also
1098 --------
1099 hermval, hermval2d, hermgrid2d, hermgrid3d
1100
1101 Notes
1102 -----
1103
1104 .. versionadded::1.7.0
1105
1106 """
1107 try:
1108 x, y, z = np.array((x, y, z), copy=0)
1109 except:
1110 raise ValueError('x, y, z are incompatible')
1111
1112 c = hermval(x, c)
1113 c = hermval(y, c, tensor=False)
1114 c = hermval(z, c, tensor=False)
1115 return c
1116
1117
1118 def hermgrid3d(x, y, z, c):
1119 """
1120 Evaluate a 3-D Hermite series on the Cartesian product of x, y, and z.
1121
1122 This function returns the values:
1123
1124 .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * H_i(a) * H_j(b) * H_k(c)
1125
1126 where the points `(a, b, c)` consist of all triples formed by taking
1127 `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
1128 a grid with `x` in the first dimension, `y` in the second, and `z` in
1129 the third.
1130
1131 The parameters `x`, `y`, and `z` are converted to arrays only if they
1132 are tuples or a lists, otherwise they are treated as a scalars. In
1133 either case, either `x`, `y`, and `z` or their elements must support
1134 multiplication and addition both with themselves and with the elements
1135 of `c`.
1136
1137 If `c` has fewer than three dimensions, ones are implicitly appended to
1138 its shape to make it 3-D. The shape of the result will be c.shape[3:] +
1139 x.shape + y.shape + z.shape.
1140
1141 Parameters
1142 ----------
1143 x, y, z : array_like, compatible objects
1144 The three dimensional series is evaluated at the points in the
1145 Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a
1146 list or tuple, it is first converted to an ndarray, otherwise it is
1147 left unchanged and, if it isn't an ndarray, it is treated as a
1148 scalar.
1149 c : array_like
1150 Array of coefficients ordered so that the coefficients for terms of
1151 degree i,j are contained in ``c[i,j]``. If `c` has dimension
1152 greater than two the remaining indices enumerate multiple sets of
1153 coefficients.
1154
1155 Returns
1156 -------
1157 values : ndarray, compatible object
1158 The values of the two dimensional polynomial at points in the Cartesian
1159 product of `x` and `y`.
1160
1161 See Also
1162 --------
1163 hermval, hermval2d, hermgrid2d, hermval3d
1164
1165 Notes
1166 -----
1167
1168 .. versionadded::1.7.0
1169
1170 """
1171 c = hermval(x, c)
1172 c = hermval(y, c)
1173 c = hermval(z, c)
1174 return c
1175
1176
1177 def hermvander(x, deg):
1178 """Pseudo-Vandermonde matrix of given degree.
1179
1180 Returns the pseudo-Vandermonde matrix of degree `deg` and sample points
1181 `x`. The pseudo-Vandermonde matrix is defined by
1182
1183 .. math:: V[..., i] = H_i(x),
1184
1185 where `0 <= i <= deg`. The leading indices of `V` index the elements of
1186 `x` and the last index is the degree of the Hermite polynomial.
1187
1188 If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the
1189 array ``V = hermvander(x, n)``, then ``np.dot(V, c)`` and
1190 ``hermval(x, c)`` are the same up to roundoff. This equivalence is
1191 useful both for least squares fitting and for the evaluation of a large
1192 number of Hermite series of the same degree and sample points.
1193
1194 Parameters
1195 ----------
1196 x : array_like
1197 Array of points. The dtype is converted to float64 or complex128
1198 depending on whether any of the elements are complex. If `x` is
1199 scalar it is converted to a 1-D array.
1200 deg : int
1201 Degree of the resulting matrix.
1202
1203 Returns
1204 -------
1205 vander : ndarray
1206 The pseudo-Vandermonde matrix. The shape of the returned matrix is
1207 ``x.shape + (deg + 1,)``, where The last index is the degree of the
1208 corresponding Hermite polynomial. The dtype will be the same as
1209 the converted `x`.
1210
1211 Examples
1212 --------
1213 >>> from numpy.polynomial.hermite import hermvander
1214 >>> x = np.array([-1, 0, 1])
1215 >>> hermvander(x, 3)
1216 array([[ 1., -2., 2., 4.],
1217 [ 1., 0., -2., -0.],
1218 [ 1., 2., 2., -4.]])
1219
1220 """
1221 ideg = int(deg)
1222 if ideg != deg:
1223 raise ValueError("deg must be integer")
1224 if ideg < 0:
1225 raise ValueError("deg must be non-negative")
1226
1227 x = np.array(x, copy=0, ndmin=1) + 0.0
1228 dims = (ideg + 1,) + x.shape
1229 dtyp = x.dtype
1230 v = np.empty(dims, dtype=dtyp)
1231 v[0] = x*0 + 1
1232 if ideg > 0:
1233 x2 = x*2
1234 v[1] = x2
1235 for i in range(2, ideg + 1):
1236 v[i] = (v[i-1]*x2 - v[i-2]*(2*(i - 1)))
1237 return np.rollaxis(v, 0, v.ndim)
1238
1239
1240 def hermvander2d(x, y, deg):
1241 """Pseudo-Vandermonde matrix of given degrees.
1242
1243 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1244 points `(x, y)`. The pseudo-Vandermonde matrix is defined by
1245
1246 .. math:: V[..., deg[1]*i + j] = H_i(x) * H_j(y),
1247
1248 where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of
1249 `V` index the points `(x, y)` and the last index encodes the degrees of
1250 the Hermite polynomials.
1251
1252 If ``V = hermvander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
1253 correspond to the elements of a 2-D coefficient array `c` of shape
1254 (xdeg + 1, ydeg + 1) in the order
1255
1256 .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ...
1257
1258 and ``np.dot(V, c.flat)`` and ``hermval2d(x, y, c)`` will be the same
1259 up to roundoff. This equivalence is useful both for least squares
1260 fitting and for the evaluation of a large number of 2-D Hermite
1261 series of the same degrees and sample points.
1262
1263 Parameters
1264 ----------
1265 x, y : array_like
1266 Arrays of point coordinates, all of the same shape. The dtypes
1267 will be converted to either float64 or complex128 depending on
1268 whether any of the elements are complex. Scalars are converted to 1-D
1269 arrays.
1270 deg : list of ints
1271 List of maximum degrees of the form [x_deg, y_deg].
1272
1273 Returns
1274 -------
1275 vander2d : ndarray
1276 The shape of the returned matrix is ``x.shape + (order,)``, where
1277 :math:`order = (deg[0]+1)*(deg([1]+1)`. The dtype will be the same
1278 as the converted `x` and `y`.
1279
1280 See Also
1281 --------
1282 hermvander, hermvander3d. hermval2d, hermval3d
1283
1284 Notes
1285 -----
1286
1287 .. versionadded::1.7.0
1288
1289 """
1290 ideg = [int(d) for d in deg]
1291 is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
1292 if is_valid != [1, 1]:
1293 raise ValueError("degrees must be non-negative integers")
1294 degx, degy = ideg
1295 x, y = np.array((x, y), copy=0) + 0.0
1296
1297 vx = hermvander(x, degx)
1298 vy = hermvander(y, degy)
1299 v = vx[..., None]*vy[..., None,:]
1300 return v.reshape(v.shape[:-2] + (-1,))
1301
1302
1303 def hermvander3d(x, y, z, deg):
1304 """Pseudo-Vandermonde matrix of given degrees.
1305
1306 Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
1307 points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`,
1308 then The pseudo-Vandermonde matrix is defined by
1309
1310 .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = H_i(x)*H_j(y)*H_k(z),
1311
1312 where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading
1313 indices of `V` index the points `(x, y, z)` and the last index encodes
1314 the degrees of the Hermite polynomials.
1315
1316 If ``V = hermvander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
1317 of `V` correspond to the elements of a 3-D coefficient array `c` of
1318 shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
1319
1320 .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
1321
1322 and ``np.dot(V, c.flat)`` and ``hermval3d(x, y, z, c)`` will be the
1323 same up to roundoff. This equivalence is useful both for least squares
1324 fitting and for the evaluation of a large number of 3-D Hermite
1325 series of the same degrees and sample points.
1326
1327 Parameters
1328 ----------
1329 x, y, z : array_like
1330 Arrays of point coordinates, all of the same shape. The dtypes will
1331 be converted to either float64 or complex128 depending on whether
1332 any of the elements are complex. Scalars are converted to 1-D
1333 arrays.
1334 deg : list of ints
1335 List of maximum degrees of the form [x_deg, y_deg, z_deg].
1336
1337 Returns
1338 -------
1339 vander3d : ndarray
1340 The shape of the returned matrix is ``x.shape + (order,)``, where
1341 :math:`order = (deg[0]+1)*(deg([1]+1)*(deg[2]+1)`. The dtype will
1342 be the same as the converted `x`, `y`, and `z`.
1343
1344 See Also
1345 --------
1346 hermvander, hermvander3d. hermval2d, hermval3d
1347
1348 Notes
1349 -----
1350
1351 .. versionadded::1.7.0
1352
1353 """
1354 ideg = [int(d) for d in deg]
1355 is_valid = [id == d and id >= 0 for id, d in zip(ideg, deg)]
1356 if is_valid != [1, 1, 1]:
1357 raise ValueError("degrees must be non-negative integers")
1358 degx, degy, degz = ideg
1359 x, y, z = np.array((x, y, z), copy=0) + 0.0
1360
1361 vx = hermvander(x, degx)
1362 vy = hermvander(y, degy)
1363 vz = hermvander(z, degz)
1364 v = vx[..., None, None]*vy[..., None,:, None]*vz[..., None, None,:]
1365 return v.reshape(v.shape[:-3] + (-1,))
1366
1367
1368 def hermfit(x, y, deg, rcond=None, full=False, w=None):
1369 """
1370 Least squares fit of Hermite series to data.
1371
1372 Return the coefficients of a Hermite series of degree `deg` that is the
1373 least squares fit to the data values `y` given at points `x`. If `y` is
1374 1-D the returned coefficients will also be 1-D. If `y` is 2-D multiple
1375 fits are done, one for each column of `y`, and the resulting
1376 coefficients are stored in the corresponding columns of a 2-D return.
1377 The fitted polynomial(s) are in the form
1378
1379 .. math:: p(x) = c_0 + c_1 * H_1(x) + ... + c_n * H_n(x),
1380
1381 where `n` is `deg`.
1382
1383 Parameters
1384 ----------
1385 x : array_like, shape (M,)
1386 x-coordinates of the M sample points ``(x[i], y[i])``.
1387 y : array_like, shape (M,) or (M, K)
1388 y-coordinates of the sample points. Several data sets of sample
1389 points sharing the same x-coordinates can be fitted at once by
1390 passing in a 2D-array that contains one dataset per column.
1391 deg : int
1392 Degree of the fitting polynomial
1393 rcond : float, optional
1394 Relative condition number of the fit. Singular values smaller than
1395 this relative to the largest singular value will be ignored. The
1396 default value is len(x)*eps, where eps is the relative precision of
1397 the float type, about 2e-16 in most cases.
1398 full : bool, optional
1399 Switch determining nature of return value. When it is False (the
1400 default) just the coefficients are returned, when True diagnostic
1401 information from the singular value decomposition is also returned.
1402 w : array_like, shape (`M`,), optional
1403 Weights. If not None, the contribution of each point
1404 ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the
1405 weights are chosen so that the errors of the products ``w[i]*y[i]``
1406 all have the same variance. The default value is None.
1407
1408 Returns
1409 -------
1410 coef : ndarray, shape (M,) or (M, K)
1411 Hermite coefficients ordered from low to high. If `y` was 2-D,
1412 the coefficients for the data in column k of `y` are in column
1413 `k`.
1414
1415 [residuals, rank, singular_values, rcond] : list
1416 These values are only returned if `full` = True
1417
1418 resid -- sum of squared residuals of the least squares fit
1419 rank -- the numerical rank of the scaled Vandermonde matrix
1420 sv -- singular values of the scaled Vandermonde matrix
1421 rcond -- value of `rcond`.
1422
1423 For more details, see `linalg.lstsq`.
1424
1425 Warns
1426 -----
1427 RankWarning
1428 The rank of the coefficient matrix in the least-squares fit is
1429 deficient. The warning is only raised if `full` = False. The
1430 warnings can be turned off by
1431
1432 >>> import warnings
1433 >>> warnings.simplefilter('ignore', RankWarning)
1434
1435 See Also
1436 --------
1437 chebfit, legfit, lagfit, polyfit, hermefit
1438 hermval : Evaluates a Hermite series.
1439 hermvander : Vandermonde matrix of Hermite series.
1440 hermweight : Hermite weight function
1441 linalg.lstsq : Computes a least-squares fit from the matrix.
1442 scipy.interpolate.UnivariateSpline : Computes spline fits.
1443
1444 Notes
1445 -----
1446 The solution is the coefficients of the Hermite series `p` that
1447 minimizes the sum of the weighted squared errors
1448
1449 .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
1450
1451 where the :math:`w_j` are the weights. This problem is solved by
1452 setting up the (typically) overdetermined matrix equation
1453
1454 .. math:: V(x) * c = w * y,
1455
1456 where `V` is the weighted pseudo Vandermonde matrix of `x`, `c` are the
1457 coefficients to be solved for, `w` are the weights, `y` are the
1458 observed values. This equation is then solved using the singular value
1459 decomposition of `V`.
1460
1461 If some of the singular values of `V` are so small that they are
1462 neglected, then a `RankWarning` will be issued. This means that the
1463 coefficient values may be poorly determined. Using a lower order fit
1464 will usually get rid of the warning. The `rcond` parameter can also be
1465 set to a value smaller than its default, but the resulting fit may be
1466 spurious and have large contributions from roundoff error.
1467
1468 Fits using Hermite series are probably most useful when the data can be
1469 approximated by ``sqrt(w(x)) * p(x)``, where `w(x)` is the Hermite
1470 weight. In that case the weight ``sqrt(w(x[i])`` should be used
1471 together with data values ``y[i]/sqrt(w(x[i])``. The weight function is
1472 available as `hermweight`.
1473
1474 References
1475 ----------
1476 .. [1] Wikipedia, "Curve fitting",
1477 http://en.wikipedia.org/wiki/Curve_fitting
1478
1479 Examples
1480 --------
1481 >>> from numpy.polynomial.hermite import hermfit, hermval
1482 >>> x = np.linspace(-10, 10)
1483 >>> err = np.random.randn(len(x))/10
1484 >>> y = hermval(x, [1, 2, 3]) + err
1485 >>> hermfit(x, y, 2)
1486 array([ 0.97902637, 1.99849131, 3.00006 ])
1487
1488 """
1489 order = int(deg) + 1
1490 x = np.asarray(x) + 0.0
1491 y = np.asarray(y) + 0.0
1492
1493 # check arguments.
1494 if deg < 0:
1495 raise ValueError("expected deg >= 0")
1496 if x.ndim != 1:
1497 raise TypeError("expected 1D vector for x")
1498 if x.size == 0:
1499 raise TypeError("expected non-empty vector for x")
1500 if y.ndim < 1 or y.ndim > 2:
1501 raise TypeError("expected 1D or 2D array for y")
1502 if len(x) != len(y):
1503 raise TypeError("expected x and y to have same length")
1504
1505 # set up the least squares matrices in transposed form
1506 lhs = hermvander(x, deg).T
1507 rhs = y.T
1508 if w is not None:
1509 w = np.asarray(w) + 0.0
1510 if w.ndim != 1:
1511 raise TypeError("expected 1D vector for w")
1512 if len(x) != len(w):
1513 raise TypeError("expected x and w to have same length")
1514 # apply weights. Don't use inplace operations as they
1515 # can cause problems with NA.
1516 lhs = lhs * w
1517 rhs = rhs * w
1518
1519 # set rcond
1520 if rcond is None:
1521 rcond = len(x)*np.finfo(x.dtype).eps
1522
1523 # Determine the norms of the design matrix columns.
1524 if issubclass(lhs.dtype.type, np.complexfloating):
1525 scl = np.sqrt((np.square(lhs.real) + np.square(lhs.imag)).sum(1))
1526 else:
1527 scl = np.sqrt(np.square(lhs).sum(1))
1528 scl[scl == 0] = 1
1529
1530 # Solve the least squares problem.
1531 c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
1532 c = (c.T/scl).T
1533
1534 # warn on rank reduction
1535 if rank != order and not full:
1536 msg = "The fit may be poorly conditioned"
1537 warnings.warn(msg, pu.RankWarning)
1538
1539 if full:
1540 return c, [resids, rank, s, rcond]
1541 else:
1542 return c
1543
1544
1545 def hermcompanion(c):
1546 """Return the scaled companion matrix of c.
1547
1548 The basis polynomials are scaled so that the companion matrix is
1549 symmetric when `c` is an Hermite basis polynomial. This provides
1550 better eigenvalue estimates than the unscaled case and for basis
1551 polynomials the eigenvalues are guaranteed to be real if
1552 `numpy.linalg.eigvalsh` is used to obtain them.
1553
1554 Parameters
1555 ----------
1556 c : array_like
1557 1-D array of Hermite series coefficients ordered from low to high
1558 degree.
1559
1560 Returns
1561 -------
1562 mat : ndarray
1563 Scaled companion matrix of dimensions (deg, deg).
1564
1565 Notes
1566 -----
1567
1568 .. versionadded::1.7.0
1569
1570 """
1571 # c is a trimmed copy
1572 [c] = pu.as_series([c])
1573 if len(c) < 2:
1574 raise ValueError('Series must have maximum degree of at least 1.')
1575 if len(c) == 2:
1576 return np.array([[-.5*c[0]/c[1]]])
1577
1578 n = len(c) - 1
1579 mat = np.zeros((n, n), dtype=c.dtype)
1580 scl = np.hstack((1., np.sqrt(2.*np.arange(1, n))))
1581 scl = np.multiply.accumulate(scl)
1582 top = mat.reshape(-1)[1::n+1]
1583 bot = mat.reshape(-1)[n::n+1]
1584 top[...] = np.sqrt(.5*np.arange(1, n))
1585 bot[...] = top
1586 mat[:, -1] -= (c[:-1]/c[-1])*(scl/scl[-1])*.5
1587 return mat
1588
1589
1590 def hermroots(c):
1591 """
1592 Compute the roots of a Hermite series.
1593
1594 Return the roots (a.k.a. "zeros") of the polynomial
1595
1596 .. math:: p(x) = \\sum_i c[i] * H_i(x).
1597
1598 Parameters
1599 ----------
1600 c : 1-D array_like
1601 1-D array of coefficients.
1602
1603 Returns
1604 -------
1605 out : ndarray
1606 Array of the roots of the series. If all the roots are real,
1607 then `out` is also real, otherwise it is complex.
1608
1609 See Also
1610 --------
1611 polyroots, legroots, lagroots, chebroots, hermeroots
1612
1613 Notes
1614 -----
1615 The root estimates are obtained as the eigenvalues of the companion
1616 matrix, Roots far from the origin of the complex plane may have large
1617 errors due to the numerical instability of the series for such
1618 values. Roots with multiplicity greater than 1 will also show larger
1619 errors as the value of the series near such points is relatively
1620 insensitive to errors in the roots. Isolated roots near the origin can
1621 be improved by a few iterations of Newton's method.
1622
1623 The Hermite series basis polynomials aren't powers of `x` so the
1624 results of this function may seem unintuitive.
1625
1626 Examples
1627 --------
1628 >>> from numpy.polynomial.hermite import hermroots, hermfromroots
1629 >>> coef = hermfromroots([-1, 0, 1])
1630 >>> coef
1631 array([ 0. , 0.25 , 0. , 0.125])
1632 >>> hermroots(coef)
1633 array([ -1.00000000e+00, -1.38777878e-17, 1.00000000e+00])
1634
1635 """
1636 # c is a trimmed copy
1637 [c] = pu.as_series([c])
1638 if len(c) <= 1:
1639 return np.array([], dtype=c.dtype)
1640 if len(c) == 2:
1641 return np.array([-.5*c[0]/c[1]])
1642
1643 m = hermcompanion(c)
1644 r = la.eigvals(m)
1645 r.sort()
1646 return r
1647
1648
1649 def hermgauss(deg):
1650 """
1651 Gauss-Hermite quadrature.
1652
1653 Computes the sample points and weights for Gauss-Hermite quadrature.
1654 These sample points and weights will correctly integrate polynomials of
1655 degree :math:`2*deg - 1` or less over the interval :math:`[-\inf, \inf]`
1656 with the weight function :math:`f(x) = \exp(-x^2)`.
1657
1658 Parameters
1659 ----------
1660 deg : int
1661 Number of sample points and weights. It must be >= 1.
1662
1663 Returns
1664 -------
1665 x : ndarray
1666 1-D ndarray containing the sample points.
1667 y : ndarray
1668 1-D ndarray containing the weights.
1669
1670 Notes
1671 -----
1672
1673 .. versionadded::1.7.0
1674
1675 The results have only been tested up to degree 100, higher degrees may
1676 be problematic. The weights are determined by using the fact that
1677
1678 .. math:: w_k = c / (H'_n(x_k) * H_{n-1}(x_k))
1679
1680 where :math:`c` is a constant independent of :math:`k` and :math:`x_k`
1681 is the k'th root of :math:`H_n`, and then scaling the results to get
1682 the right value when integrating 1.
1683
1684 """
1685 ideg = int(deg)
1686 if ideg != deg or ideg < 1:
1687 raise ValueError("deg must be a non-negative integer")
1688
1689 # first approximation of roots. We use the fact that the companion
1690 # matrix is symmetric in this case in order to obtain better zeros.
1691 c = np.array([0]*deg + [1])
1692 m = hermcompanion(c)
1693 x = la.eigvals(m)
1694 x.sort()
1695
1696 # improve roots by one application of Newton
1697 dy = hermval(x, c)
1698 df = hermval(x, hermder(c))
1699 x -= dy/df
1700
1701 # compute the weights. We scale the factor to avoid possible numerical
1702 # overflow.
1703 fm = hermval(x, c[1:])
1704 fm /= np.abs(fm).max()
1705 df /= np.abs(df).max()
1706 w = 1/(fm * df)
1707
1708 # for Hermite we can also symmetrize
1709 w = (w + w[::-1])/2
1710 x = (x - x[::-1])/2
1711
1712 # scale w to get the right value
1713 w *= np.sqrt(np.pi) / w.sum()
1714
1715 return x, w
1716
1717
1718 def hermweight(x):
1719 """
1720 Weight function of the Hermite polynomials.
1721
1722 The weight function is :math:`\exp(-x^2)` and the interval of
1723 integration is :math:`[-\inf, \inf]`. the Hermite polynomials are
1724 orthogonal, but not normalized, with respect to this weight function.
1725
1726 Parameters
1727 ----------
1728 x : array_like
1729 Values at which the weight function will be computed.
1730
1731 Returns
1732 -------
1733 w : ndarray
1734 The weight function at `x`.
1735
1736 Notes
1737 -----
1738
1739 .. versionadded::1.7.0
1740
1741 """
1742 w = np.exp(-x**2)
1743 return w
1744
1745
1746 #
1747 # Hermite series class
1748 #
1749
1750 class Hermite(ABCPolyBase):
1751 """An Hermite series class.
1752
1753 The Hermite class provides the standard Python numerical methods
1754 '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
1755 attributes and methods listed in the `ABCPolyBase` documentation.
1756
1757 Parameters
1758 ----------
1759 coef : array_like
1760 Laguerre coefficients in order of increasing degree, i.e,
1761 ``(1, 2, 3)`` gives ``1*H_0(x) + 2*H_1(X) + 3*H_2(x)``.
1762 domain : (2,) array_like, optional
1763 Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
1764 to the interval ``[window[0], window[1]]`` by shifting and scaling.
1765 The default value is [-1, 1].
1766 window : (2,) array_like, optional
1767 Window, see `domain` for its use. The default value is [-1, 1].
1768
1769 .. versionadded:: 1.6.0
1770
1771 """
1772 # Virtual Functions
1773 _add = staticmethod(hermadd)
1774 _sub = staticmethod(hermsub)
1775 _mul = staticmethod(hermmul)
1776 _div = staticmethod(hermdiv)
1777 _pow = staticmethod(hermpow)
1778 _val = staticmethod(hermval)
1779 _int = staticmethod(hermint)
1780 _der = staticmethod(hermder)
1781 _fit = staticmethod(hermfit)
1782 _line = staticmethod(hermline)
1783 _roots = staticmethod(hermroots)
1784 _fromroots = staticmethod(hermfromroots)
1785
1786 # Virtual properties
1787 nickname = 'herm'
1788 domain = np.array(hermdomain)
1789 window = np.array(hermdomain)