hermite_e.py 50 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675
  1. """
  2. ===================================================================
  3. HermiteE Series, "Probabilists" (:mod:`numpy.polynomial.hermite_e`)
  4. ===================================================================
  5. This module provides a number of objects (mostly functions) useful for
  6. dealing with Hermite_e series, including a `HermiteE` class that
  7. encapsulates the usual arithmetic operations. (General information
  8. on how this module represents and works with such polynomials is in the
  9. docstring for its "parent" sub-package, `numpy.polynomial`).
  10. Classes
  11. -------
  12. .. autosummary::
  13. :toctree: generated/
  14. HermiteE
  15. Constants
  16. ---------
  17. .. autosummary::
  18. :toctree: generated/
  19. hermedomain
  20. hermezero
  21. hermeone
  22. hermex
  23. Arithmetic
  24. ----------
  25. .. autosummary::
  26. :toctree: generated/
  27. hermeadd
  28. hermesub
  29. hermemulx
  30. hermemul
  31. hermediv
  32. hermepow
  33. hermeval
  34. hermeval2d
  35. hermeval3d
  36. hermegrid2d
  37. hermegrid3d
  38. Calculus
  39. --------
  40. .. autosummary::
  41. :toctree: generated/
  42. hermeder
  43. hermeint
  44. Misc Functions
  45. --------------
  46. .. autosummary::
  47. :toctree: generated/
  48. hermefromroots
  49. hermeroots
  50. hermevander
  51. hermevander2d
  52. hermevander3d
  53. hermegauss
  54. hermeweight
  55. hermecompanion
  56. hermefit
  57. hermetrim
  58. hermeline
  59. herme2poly
  60. poly2herme
  61. See also
  62. --------
  63. `numpy.polynomial`
  64. """
  65. import numpy as np
  66. import numpy.linalg as la
  67. from numpy.core.multiarray import normalize_axis_index
  68. from . import polyutils as pu
  69. from ._polybase import ABCPolyBase
  70. __all__ = [
  71. 'hermezero', 'hermeone', 'hermex', 'hermedomain', 'hermeline',
  72. 'hermeadd', 'hermesub', 'hermemulx', 'hermemul', 'hermediv',
  73. 'hermepow', 'hermeval', 'hermeder', 'hermeint', 'herme2poly',
  74. 'poly2herme', 'hermefromroots', 'hermevander', 'hermefit', 'hermetrim',
  75. 'hermeroots', 'HermiteE', 'hermeval2d', 'hermeval3d', 'hermegrid2d',
  76. 'hermegrid3d', 'hermevander2d', 'hermevander3d', 'hermecompanion',
  77. 'hermegauss', 'hermeweight']
  78. hermetrim = pu.trimcoef
  79. def poly2herme(pol):
  80. """
  81. poly2herme(pol)
  82. Convert a polynomial to a Hermite series.
  83. Convert an array representing the coefficients of a polynomial (relative
  84. to the "standard" basis) ordered from lowest degree to highest, to an
  85. array of the coefficients of the equivalent Hermite series, ordered
  86. from lowest to highest degree.
  87. Parameters
  88. ----------
  89. pol : array_like
  90. 1-D array containing the polynomial coefficients
  91. Returns
  92. -------
  93. c : ndarray
  94. 1-D array containing the coefficients of the equivalent Hermite
  95. series.
  96. See Also
  97. --------
  98. herme2poly
  99. Notes
  100. -----
  101. The easy way to do conversions between polynomial basis sets
  102. is to use the convert method of a class instance.
  103. Examples
  104. --------
  105. >>> from numpy.polynomial.hermite_e import poly2herme
  106. >>> poly2herme(np.arange(4))
  107. array([ 2., 10., 2., 3.])
  108. """
  109. [pol] = pu.as_series([pol])
  110. deg = len(pol) - 1
  111. res = 0
  112. for i in range(deg, -1, -1):
  113. res = hermeadd(hermemulx(res), pol[i])
  114. return res
  115. def herme2poly(c):
  116. """
  117. Convert a Hermite series to a polynomial.
  118. Convert an array representing the coefficients of a Hermite series,
  119. ordered from lowest degree to highest, to an array of the coefficients
  120. of the equivalent polynomial (relative to the "standard" basis) ordered
  121. from lowest to highest degree.
  122. Parameters
  123. ----------
  124. c : array_like
  125. 1-D array containing the Hermite series coefficients, ordered
  126. from lowest order term to highest.
  127. Returns
  128. -------
  129. pol : ndarray
  130. 1-D array containing the coefficients of the equivalent polynomial
  131. (relative to the "standard" basis) ordered from lowest order term
  132. to highest.
  133. See Also
  134. --------
  135. poly2herme
  136. Notes
  137. -----
  138. The easy way to do conversions between polynomial basis sets
  139. is to use the convert method of a class instance.
  140. Examples
  141. --------
  142. >>> from numpy.polynomial.hermite_e import herme2poly
  143. >>> herme2poly([ 2., 10., 2., 3.])
  144. array([0., 1., 2., 3.])
  145. """
  146. from .polynomial import polyadd, polysub, polymulx
  147. [c] = pu.as_series([c])
  148. n = len(c)
  149. if n == 1:
  150. return c
  151. if n == 2:
  152. return c
  153. else:
  154. c0 = c[-2]
  155. c1 = c[-1]
  156. # i is the current degree of c1
  157. for i in range(n - 1, 1, -1):
  158. tmp = c0
  159. c0 = polysub(c[i - 2], c1*(i - 1))
  160. c1 = polyadd(tmp, polymulx(c1))
  161. return polyadd(c0, polymulx(c1))
  162. #
  163. # These are constant arrays are of integer type so as to be compatible
  164. # with the widest range of other types, such as Decimal.
  165. #
  166. # Hermite
  167. hermedomain = np.array([-1, 1])
  168. # Hermite coefficients representing zero.
  169. hermezero = np.array([0])
  170. # Hermite coefficients representing one.
  171. hermeone = np.array([1])
  172. # Hermite coefficients representing the identity x.
  173. hermex = np.array([0, 1])
  174. def hermeline(off, scl):
  175. """
  176. Hermite series whose graph is a straight line.
  177. Parameters
  178. ----------
  179. off, scl : scalars
  180. The specified line is given by ``off + scl*x``.
  181. Returns
  182. -------
  183. y : ndarray
  184. This module's representation of the Hermite series for
  185. ``off + scl*x``.
  186. See Also
  187. --------
  188. polyline, chebline
  189. Examples
  190. --------
  191. >>> from numpy.polynomial.hermite_e import hermeline
  192. >>> from numpy.polynomial.hermite_e import hermeline, hermeval
  193. >>> hermeval(0,hermeline(3, 2))
  194. 3.0
  195. >>> hermeval(1,hermeline(3, 2))
  196. 5.0
  197. """
  198. if scl != 0:
  199. return np.array([off, scl])
  200. else:
  201. return np.array([off])
  202. def hermefromroots(roots):
  203. """
  204. Generate a HermiteE series with given roots.
  205. The function returns the coefficients of the polynomial
  206. .. math:: p(x) = (x - r_0) * (x - r_1) * ... * (x - r_n),
  207. in HermiteE form, where the `r_n` are the roots specified in `roots`.
  208. If a zero has multiplicity n, then it must appear in `roots` n times.
  209. For instance, if 2 is a root of multiplicity three and 3 is a root of
  210. multiplicity 2, then `roots` looks something like [2, 2, 2, 3, 3]. The
  211. roots can appear in any order.
  212. If the returned coefficients are `c`, then
  213. .. math:: p(x) = c_0 + c_1 * He_1(x) + ... + c_n * He_n(x)
  214. The coefficient of the last term is not generally 1 for monic
  215. polynomials in HermiteE form.
  216. Parameters
  217. ----------
  218. roots : array_like
  219. Sequence containing the roots.
  220. Returns
  221. -------
  222. out : ndarray
  223. 1-D array of coefficients. If all roots are real then `out` is a
  224. real array, if some of the roots are complex, then `out` is complex
  225. even if all the coefficients in the result are real (see Examples
  226. below).
  227. See Also
  228. --------
  229. polyfromroots, legfromroots, lagfromroots, hermfromroots, chebfromroots
  230. Examples
  231. --------
  232. >>> from numpy.polynomial.hermite_e import hermefromroots, hermeval
  233. >>> coef = hermefromroots((-1, 0, 1))
  234. >>> hermeval((-1, 0, 1), coef)
  235. array([0., 0., 0.])
  236. >>> coef = hermefromroots((-1j, 1j))
  237. >>> hermeval((-1j, 1j), coef)
  238. array([0.+0.j, 0.+0.j])
  239. """
  240. return pu._fromroots(hermeline, hermemul, roots)
  241. def hermeadd(c1, c2):
  242. """
  243. Add one Hermite series to another.
  244. Returns the sum of two Hermite series `c1` + `c2`. The arguments
  245. are sequences of coefficients ordered from lowest order term to
  246. highest, i.e., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
  247. Parameters
  248. ----------
  249. c1, c2 : array_like
  250. 1-D arrays of Hermite series coefficients ordered from low to
  251. high.
  252. Returns
  253. -------
  254. out : ndarray
  255. Array representing the Hermite series of their sum.
  256. See Also
  257. --------
  258. hermesub, hermemulx, hermemul, hermediv, hermepow
  259. Notes
  260. -----
  261. Unlike multiplication, division, etc., the sum of two Hermite series
  262. is a Hermite series (without having to "reproject" the result onto
  263. the basis set) so addition, just like that of "standard" polynomials,
  264. is simply "component-wise."
  265. Examples
  266. --------
  267. >>> from numpy.polynomial.hermite_e import hermeadd
  268. >>> hermeadd([1, 2, 3], [1, 2, 3, 4])
  269. array([2., 4., 6., 4.])
  270. """
  271. return pu._add(c1, c2)
  272. def hermesub(c1, c2):
  273. """
  274. Subtract one Hermite series from another.
  275. Returns the difference of two Hermite series `c1` - `c2`. The
  276. sequences of coefficients are from lowest order term to highest, i.e.,
  277. [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
  278. Parameters
  279. ----------
  280. c1, c2 : array_like
  281. 1-D arrays of Hermite series coefficients ordered from low to
  282. high.
  283. Returns
  284. -------
  285. out : ndarray
  286. Of Hermite series coefficients representing their difference.
  287. See Also
  288. --------
  289. hermeadd, hermemulx, hermemul, hermediv, hermepow
  290. Notes
  291. -----
  292. Unlike multiplication, division, etc., the difference of two Hermite
  293. series is a Hermite series (without having to "reproject" the result
  294. onto the basis set) so subtraction, just like that of "standard"
  295. polynomials, is simply "component-wise."
  296. Examples
  297. --------
  298. >>> from numpy.polynomial.hermite_e import hermesub
  299. >>> hermesub([1, 2, 3, 4], [1, 2, 3])
  300. array([0., 0., 0., 4.])
  301. """
  302. return pu._sub(c1, c2)
  303. def hermemulx(c):
  304. """Multiply a Hermite series by x.
  305. Multiply the Hermite series `c` by x, where x is the independent
  306. variable.
  307. Parameters
  308. ----------
  309. c : array_like
  310. 1-D array of Hermite series coefficients ordered from low to
  311. high.
  312. Returns
  313. -------
  314. out : ndarray
  315. Array representing the result of the multiplication.
  316. Notes
  317. -----
  318. The multiplication uses the recursion relationship for Hermite
  319. polynomials in the form
  320. .. math::
  321. xP_i(x) = (P_{i + 1}(x) + iP_{i - 1}(x)))
  322. Examples
  323. --------
  324. >>> from numpy.polynomial.hermite_e import hermemulx
  325. >>> hermemulx([1, 2, 3])
  326. array([2., 7., 2., 3.])
  327. """
  328. # c is a trimmed copy
  329. [c] = pu.as_series([c])
  330. # The zero series needs special treatment
  331. if len(c) == 1 and c[0] == 0:
  332. return c
  333. prd = np.empty(len(c) + 1, dtype=c.dtype)
  334. prd[0] = c[0]*0
  335. prd[1] = c[0]
  336. for i in range(1, len(c)):
  337. prd[i + 1] = c[i]
  338. prd[i - 1] += c[i]*i
  339. return prd
  340. def hermemul(c1, c2):
  341. """
  342. Multiply one Hermite series by another.
  343. Returns the product of two Hermite series `c1` * `c2`. The arguments
  344. are sequences of coefficients, from lowest order "term" to highest,
  345. e.g., [1,2,3] represents the series ``P_0 + 2*P_1 + 3*P_2``.
  346. Parameters
  347. ----------
  348. c1, c2 : array_like
  349. 1-D arrays of Hermite series coefficients ordered from low to
  350. high.
  351. Returns
  352. -------
  353. out : ndarray
  354. Of Hermite series coefficients representing their product.
  355. See Also
  356. --------
  357. hermeadd, hermesub, hermemulx, hermediv, hermepow
  358. Notes
  359. -----
  360. In general, the (polynomial) product of two C-series results in terms
  361. that are not in the Hermite polynomial basis set. Thus, to express
  362. the product as a Hermite series, it is necessary to "reproject" the
  363. product onto said basis set, which may produce "unintuitive" (but
  364. correct) results; see Examples section below.
  365. Examples
  366. --------
  367. >>> from numpy.polynomial.hermite_e import hermemul
  368. >>> hermemul([1, 2, 3], [0, 1, 2])
  369. array([14., 15., 28., 7., 6.])
  370. """
  371. # s1, s2 are trimmed copies
  372. [c1, c2] = pu.as_series([c1, c2])
  373. if len(c1) > len(c2):
  374. c = c2
  375. xs = c1
  376. else:
  377. c = c1
  378. xs = c2
  379. if len(c) == 1:
  380. c0 = c[0]*xs
  381. c1 = 0
  382. elif len(c) == 2:
  383. c0 = c[0]*xs
  384. c1 = c[1]*xs
  385. else:
  386. nd = len(c)
  387. c0 = c[-2]*xs
  388. c1 = c[-1]*xs
  389. for i in range(3, len(c) + 1):
  390. tmp = c0
  391. nd = nd - 1
  392. c0 = hermesub(c[-i]*xs, c1*(nd - 1))
  393. c1 = hermeadd(tmp, hermemulx(c1))
  394. return hermeadd(c0, hermemulx(c1))
  395. def hermediv(c1, c2):
  396. """
  397. Divide one Hermite series by another.
  398. Returns the quotient-with-remainder of two Hermite series
  399. `c1` / `c2`. The arguments are sequences of coefficients from lowest
  400. order "term" to highest, e.g., [1,2,3] represents the series
  401. ``P_0 + 2*P_1 + 3*P_2``.
  402. Parameters
  403. ----------
  404. c1, c2 : array_like
  405. 1-D arrays of Hermite series coefficients ordered from low to
  406. high.
  407. Returns
  408. -------
  409. [quo, rem] : ndarrays
  410. Of Hermite series coefficients representing the quotient and
  411. remainder.
  412. See Also
  413. --------
  414. hermeadd, hermesub, hermemulx, hermemul, hermepow
  415. Notes
  416. -----
  417. In general, the (polynomial) division of one Hermite series by another
  418. results in quotient and remainder terms that are not in the Hermite
  419. polynomial basis set. Thus, to express these results as a Hermite
  420. series, it is necessary to "reproject" the results onto the Hermite
  421. basis set, which may produce "unintuitive" (but correct) results; see
  422. Examples section below.
  423. Examples
  424. --------
  425. >>> from numpy.polynomial.hermite_e import hermediv
  426. >>> hermediv([ 14., 15., 28., 7., 6.], [0, 1, 2])
  427. (array([1., 2., 3.]), array([0.]))
  428. >>> hermediv([ 15., 17., 28., 7., 6.], [0, 1, 2])
  429. (array([1., 2., 3.]), array([1., 2.]))
  430. """
  431. return pu._div(hermemul, c1, c2)
  432. def hermepow(c, pow, maxpower=16):
  433. """Raise a Hermite series to a power.
  434. Returns the Hermite series `c` raised to the power `pow`. The
  435. argument `c` is a sequence of coefficients ordered from low to high.
  436. i.e., [1,2,3] is the series ``P_0 + 2*P_1 + 3*P_2.``
  437. Parameters
  438. ----------
  439. c : array_like
  440. 1-D array of Hermite series coefficients ordered from low to
  441. high.
  442. pow : integer
  443. Power to which the series will be raised
  444. maxpower : integer, optional
  445. Maximum power allowed. This is mainly to limit growth of the series
  446. to unmanageable size. Default is 16
  447. Returns
  448. -------
  449. coef : ndarray
  450. Hermite series of power.
  451. See Also
  452. --------
  453. hermeadd, hermesub, hermemulx, hermemul, hermediv
  454. Examples
  455. --------
  456. >>> from numpy.polynomial.hermite_e import hermepow
  457. >>> hermepow([1, 2, 3], 2)
  458. array([23., 28., 46., 12., 9.])
  459. """
  460. return pu._pow(hermemul, c, pow, maxpower)
  461. def hermeder(c, m=1, scl=1, axis=0):
  462. """
  463. Differentiate a Hermite_e series.
  464. Returns the series coefficients `c` differentiated `m` times along
  465. `axis`. At each iteration the result is multiplied by `scl` (the
  466. scaling factor is for use in a linear change of variable). The argument
  467. `c` is an array of coefficients from low to high degree along each
  468. axis, e.g., [1,2,3] represents the series ``1*He_0 + 2*He_1 + 3*He_2``
  469. while [[1,2],[1,2]] represents ``1*He_0(x)*He_0(y) + 1*He_1(x)*He_0(y)
  470. + 2*He_0(x)*He_1(y) + 2*He_1(x)*He_1(y)`` if axis=0 is ``x`` and axis=1
  471. is ``y``.
  472. Parameters
  473. ----------
  474. c : array_like
  475. Array of Hermite_e series coefficients. If `c` is multidimensional
  476. the different axis correspond to different variables with the
  477. degree in each axis given by the corresponding index.
  478. m : int, optional
  479. Number of derivatives taken, must be non-negative. (Default: 1)
  480. scl : scalar, optional
  481. Each differentiation is multiplied by `scl`. The end result is
  482. multiplication by ``scl**m``. This is for use in a linear change of
  483. variable. (Default: 1)
  484. axis : int, optional
  485. Axis over which the derivative is taken. (Default: 0).
  486. .. versionadded:: 1.7.0
  487. Returns
  488. -------
  489. der : ndarray
  490. Hermite series of the derivative.
  491. See Also
  492. --------
  493. hermeint
  494. Notes
  495. -----
  496. In general, the result of differentiating a Hermite series does not
  497. resemble the same operation on a power series. Thus the result of this
  498. function may be "unintuitive," albeit correct; see Examples section
  499. below.
  500. Examples
  501. --------
  502. >>> from numpy.polynomial.hermite_e import hermeder
  503. >>> hermeder([ 1., 1., 1., 1.])
  504. array([1., 2., 3.])
  505. >>> hermeder([-0.25, 1., 1./2., 1./3., 1./4 ], m=2)
  506. array([1., 2., 3.])
  507. """
  508. c = np.array(c, ndmin=1, copy=True)
  509. if c.dtype.char in '?bBhHiIlLqQpP':
  510. c = c.astype(np.double)
  511. cnt = pu._deprecate_as_int(m, "the order of derivation")
  512. iaxis = pu._deprecate_as_int(axis, "the axis")
  513. if cnt < 0:
  514. raise ValueError("The order of derivation must be non-negative")
  515. iaxis = normalize_axis_index(iaxis, c.ndim)
  516. if cnt == 0:
  517. return c
  518. c = np.moveaxis(c, iaxis, 0)
  519. n = len(c)
  520. if cnt >= n:
  521. return c[:1]*0
  522. else:
  523. for i in range(cnt):
  524. n = n - 1
  525. c *= scl
  526. der = np.empty((n,) + c.shape[1:], dtype=c.dtype)
  527. for j in range(n, 0, -1):
  528. der[j - 1] = j*c[j]
  529. c = der
  530. c = np.moveaxis(c, 0, iaxis)
  531. return c
  532. def hermeint(c, m=1, k=[], lbnd=0, scl=1, axis=0):
  533. """
  534. Integrate a Hermite_e series.
  535. Returns the Hermite_e series coefficients `c` integrated `m` times from
  536. `lbnd` along `axis`. At each iteration the resulting series is
  537. **multiplied** by `scl` and an integration constant, `k`, is added.
  538. The scaling factor is for use in a linear change of variable. ("Buyer
  539. beware": note that, depending on what one is doing, one may want `scl`
  540. to be the reciprocal of what one might expect; for more information,
  541. see the Notes section below.) The argument `c` is an array of
  542. coefficients from low to high degree along each axis, e.g., [1,2,3]
  543. represents the series ``H_0 + 2*H_1 + 3*H_2`` while [[1,2],[1,2]]
  544. represents ``1*H_0(x)*H_0(y) + 1*H_1(x)*H_0(y) + 2*H_0(x)*H_1(y) +
  545. 2*H_1(x)*H_1(y)`` if axis=0 is ``x`` and axis=1 is ``y``.
  546. Parameters
  547. ----------
  548. c : array_like
  549. Array of Hermite_e series coefficients. If c is multidimensional
  550. the different axis correspond to different variables with the
  551. degree in each axis given by the corresponding index.
  552. m : int, optional
  553. Order of integration, must be positive. (Default: 1)
  554. k : {[], list, scalar}, optional
  555. Integration constant(s). The value of the first integral at
  556. ``lbnd`` is the first value in the list, the value of the second
  557. integral at ``lbnd`` is the second value, etc. If ``k == []`` (the
  558. default), all constants are set to zero. If ``m == 1``, a single
  559. scalar can be given instead of a list.
  560. lbnd : scalar, optional
  561. The lower bound of the integral. (Default: 0)
  562. scl : scalar, optional
  563. Following each integration the result is *multiplied* by `scl`
  564. before the integration constant is added. (Default: 1)
  565. axis : int, optional
  566. Axis over which the integral is taken. (Default: 0).
  567. .. versionadded:: 1.7.0
  568. Returns
  569. -------
  570. S : ndarray
  571. Hermite_e series coefficients of the integral.
  572. Raises
  573. ------
  574. ValueError
  575. If ``m < 0``, ``len(k) > m``, ``np.ndim(lbnd) != 0``, or
  576. ``np.ndim(scl) != 0``.
  577. See Also
  578. --------
  579. hermeder
  580. Notes
  581. -----
  582. Note that the result of each integration is *multiplied* by `scl`.
  583. Why is this important to note? Say one is making a linear change of
  584. variable :math:`u = ax + b` in an integral relative to `x`. Then
  585. :math:`dx = du/a`, so one will need to set `scl` equal to
  586. :math:`1/a` - perhaps not what one would have first thought.
  587. Also note that, in general, the result of integrating a C-series needs
  588. to be "reprojected" onto the C-series basis set. Thus, typically,
  589. the result of this function is "unintuitive," albeit correct; see
  590. Examples section below.
  591. Examples
  592. --------
  593. >>> from numpy.polynomial.hermite_e import hermeint
  594. >>> hermeint([1, 2, 3]) # integrate once, value 0 at 0.
  595. array([1., 1., 1., 1.])
  596. >>> hermeint([1, 2, 3], m=2) # integrate twice, value & deriv 0 at 0
  597. array([-0.25 , 1. , 0.5 , 0.33333333, 0.25 ]) # may vary
  598. >>> hermeint([1, 2, 3], k=1) # integrate once, value 1 at 0.
  599. array([2., 1., 1., 1.])
  600. >>> hermeint([1, 2, 3], lbnd=-1) # integrate once, value 0 at -1
  601. array([-1., 1., 1., 1.])
  602. >>> hermeint([1, 2, 3], m=2, k=[1, 2], lbnd=-1)
  603. array([ 1.83333333, 0. , 0.5 , 0.33333333, 0.25 ]) # may vary
  604. """
  605. c = np.array(c, ndmin=1, copy=True)
  606. if c.dtype.char in '?bBhHiIlLqQpP':
  607. c = c.astype(np.double)
  608. if not np.iterable(k):
  609. k = [k]
  610. cnt = pu._deprecate_as_int(m, "the order of integration")
  611. iaxis = pu._deprecate_as_int(axis, "the axis")
  612. if cnt < 0:
  613. raise ValueError("The order of integration must be non-negative")
  614. if len(k) > cnt:
  615. raise ValueError("Too many integration constants")
  616. if np.ndim(lbnd) != 0:
  617. raise ValueError("lbnd must be a scalar.")
  618. if np.ndim(scl) != 0:
  619. raise ValueError("scl must be a scalar.")
  620. iaxis = normalize_axis_index(iaxis, c.ndim)
  621. if cnt == 0:
  622. return c
  623. c = np.moveaxis(c, iaxis, 0)
  624. k = list(k) + [0]*(cnt - len(k))
  625. for i in range(cnt):
  626. n = len(c)
  627. c *= scl
  628. if n == 1 and np.all(c[0] == 0):
  629. c[0] += k[i]
  630. else:
  631. tmp = np.empty((n + 1,) + c.shape[1:], dtype=c.dtype)
  632. tmp[0] = c[0]*0
  633. tmp[1] = c[0]
  634. for j in range(1, n):
  635. tmp[j + 1] = c[j]/(j + 1)
  636. tmp[0] += k[i] - hermeval(lbnd, tmp)
  637. c = tmp
  638. c = np.moveaxis(c, 0, iaxis)
  639. return c
  640. def hermeval(x, c, tensor=True):
  641. """
  642. Evaluate an HermiteE series at points x.
  643. If `c` is of length `n + 1`, this function returns the value:
  644. .. math:: p(x) = c_0 * He_0(x) + c_1 * He_1(x) + ... + c_n * He_n(x)
  645. The parameter `x` is converted to an array only if it is a tuple or a
  646. list, otherwise it is treated as a scalar. In either case, either `x`
  647. or its elements must support multiplication and addition both with
  648. themselves and with the elements of `c`.
  649. If `c` is a 1-D array, then `p(x)` will have the same shape as `x`. If
  650. `c` is multidimensional, then the shape of the result depends on the
  651. value of `tensor`. If `tensor` is true the shape will be c.shape[1:] +
  652. x.shape. If `tensor` is false the shape will be c.shape[1:]. Note that
  653. scalars have shape (,).
  654. Trailing zeros in the coefficients will be used in the evaluation, so
  655. they should be avoided if efficiency is a concern.
  656. Parameters
  657. ----------
  658. x : array_like, compatible object
  659. If `x` is a list or tuple, it is converted to an ndarray, otherwise
  660. it is left unchanged and treated as a scalar. In either case, `x`
  661. or its elements must support addition and multiplication with
  662. with themselves and with the elements of `c`.
  663. c : array_like
  664. Array of coefficients ordered so that the coefficients for terms of
  665. degree n are contained in c[n]. If `c` is multidimensional the
  666. remaining indices enumerate multiple polynomials. In the two
  667. dimensional case the coefficients may be thought of as stored in
  668. the columns of `c`.
  669. tensor : boolean, optional
  670. If True, the shape of the coefficient array is extended with ones
  671. on the right, one for each dimension of `x`. Scalars have dimension 0
  672. for this action. The result is that every column of coefficients in
  673. `c` is evaluated for every element of `x`. If False, `x` is broadcast
  674. over the columns of `c` for the evaluation. This keyword is useful
  675. when `c` is multidimensional. The default value is True.
  676. .. versionadded:: 1.7.0
  677. Returns
  678. -------
  679. values : ndarray, algebra_like
  680. The shape of the return value is described above.
  681. See Also
  682. --------
  683. hermeval2d, hermegrid2d, hermeval3d, hermegrid3d
  684. Notes
  685. -----
  686. The evaluation uses Clenshaw recursion, aka synthetic division.
  687. Examples
  688. --------
  689. >>> from numpy.polynomial.hermite_e import hermeval
  690. >>> coef = [1,2,3]
  691. >>> hermeval(1, coef)
  692. 3.0
  693. >>> hermeval([[1,2],[3,4]], coef)
  694. array([[ 3., 14.],
  695. [31., 54.]])
  696. """
  697. c = np.array(c, ndmin=1, copy=False)
  698. if c.dtype.char in '?bBhHiIlLqQpP':
  699. c = c.astype(np.double)
  700. if isinstance(x, (tuple, list)):
  701. x = np.asarray(x)
  702. if isinstance(x, np.ndarray) and tensor:
  703. c = c.reshape(c.shape + (1,)*x.ndim)
  704. if len(c) == 1:
  705. c0 = c[0]
  706. c1 = 0
  707. elif len(c) == 2:
  708. c0 = c[0]
  709. c1 = c[1]
  710. else:
  711. nd = len(c)
  712. c0 = c[-2]
  713. c1 = c[-1]
  714. for i in range(3, len(c) + 1):
  715. tmp = c0
  716. nd = nd - 1
  717. c0 = c[-i] - c1*(nd - 1)
  718. c1 = tmp + c1*x
  719. return c0 + c1*x
  720. def hermeval2d(x, y, c):
  721. """
  722. Evaluate a 2-D HermiteE series at points (x, y).
  723. This function returns the values:
  724. .. math:: p(x,y) = \\sum_{i,j} c_{i,j} * He_i(x) * He_j(y)
  725. The parameters `x` and `y` are converted to arrays only if they are
  726. tuples or a lists, otherwise they are treated as a scalars and they
  727. must have the same shape after conversion. In either case, either `x`
  728. and `y` or their elements must support multiplication and addition both
  729. with themselves and with the elements of `c`.
  730. If `c` is a 1-D array a one is implicitly appended to its shape to make
  731. it 2-D. The shape of the result will be c.shape[2:] + x.shape.
  732. Parameters
  733. ----------
  734. x, y : array_like, compatible objects
  735. The two dimensional series is evaluated at the points `(x, y)`,
  736. where `x` and `y` must have the same shape. If `x` or `y` is a list
  737. or tuple, it is first converted to an ndarray, otherwise it is left
  738. unchanged and if it isn't an ndarray it is treated as a scalar.
  739. c : array_like
  740. Array of coefficients ordered so that the coefficient of the term
  741. of multi-degree i,j is contained in ``c[i,j]``. If `c` has
  742. dimension greater than two the remaining indices enumerate multiple
  743. sets of coefficients.
  744. Returns
  745. -------
  746. values : ndarray, compatible object
  747. The values of the two dimensional polynomial at points formed with
  748. pairs of corresponding values from `x` and `y`.
  749. See Also
  750. --------
  751. hermeval, hermegrid2d, hermeval3d, hermegrid3d
  752. Notes
  753. -----
  754. .. versionadded:: 1.7.0
  755. """
  756. return pu._valnd(hermeval, c, x, y)
  757. def hermegrid2d(x, y, c):
  758. """
  759. Evaluate a 2-D HermiteE series on the Cartesian product of x and y.
  760. This function returns the values:
  761. .. math:: p(a,b) = \\sum_{i,j} c_{i,j} * H_i(a) * H_j(b)
  762. where the points `(a, b)` consist of all pairs formed by taking
  763. `a` from `x` and `b` from `y`. The resulting points form a grid with
  764. `x` in the first dimension and `y` in the second.
  765. The parameters `x` and `y` are converted to arrays only if they are
  766. tuples or a lists, otherwise they are treated as a scalars. In either
  767. case, either `x` and `y` or their elements must support multiplication
  768. and addition both with themselves and with the elements of `c`.
  769. If `c` has fewer than two dimensions, ones are implicitly appended to
  770. its shape to make it 2-D. The shape of the result will be c.shape[2:] +
  771. x.shape.
  772. Parameters
  773. ----------
  774. x, y : array_like, compatible objects
  775. The two dimensional series is evaluated at the points in the
  776. Cartesian product of `x` and `y`. If `x` or `y` is a list or
  777. tuple, it is first converted to an ndarray, otherwise it is left
  778. unchanged and, if it isn't an ndarray, it is treated as a scalar.
  779. c : array_like
  780. Array of coefficients ordered so that the coefficients for terms of
  781. degree i,j are contained in ``c[i,j]``. If `c` has dimension
  782. greater than two the remaining indices enumerate multiple sets of
  783. coefficients.
  784. Returns
  785. -------
  786. values : ndarray, compatible object
  787. The values of the two dimensional polynomial at points in the Cartesian
  788. product of `x` and `y`.
  789. See Also
  790. --------
  791. hermeval, hermeval2d, hermeval3d, hermegrid3d
  792. Notes
  793. -----
  794. .. versionadded:: 1.7.0
  795. """
  796. return pu._gridnd(hermeval, c, x, y)
  797. def hermeval3d(x, y, z, c):
  798. """
  799. Evaluate a 3-D Hermite_e series at points (x, y, z).
  800. This function returns the values:
  801. .. math:: p(x,y,z) = \\sum_{i,j,k} c_{i,j,k} * He_i(x) * He_j(y) * He_k(z)
  802. The parameters `x`, `y`, and `z` are converted to arrays only if
  803. they are tuples or a lists, otherwise they are treated as a scalars and
  804. they must have the same shape after conversion. In either case, either
  805. `x`, `y`, and `z` or their elements must support multiplication and
  806. addition both with themselves and with the elements of `c`.
  807. If `c` has fewer than 3 dimensions, ones are implicitly appended to its
  808. shape to make it 3-D. The shape of the result will be c.shape[3:] +
  809. x.shape.
  810. Parameters
  811. ----------
  812. x, y, z : array_like, compatible object
  813. The three dimensional series is evaluated at the points
  814. `(x, y, z)`, where `x`, `y`, and `z` must have the same shape. If
  815. any of `x`, `y`, or `z` is a list or tuple, it is first converted
  816. to an ndarray, otherwise it is left unchanged and if it isn't an
  817. ndarray it is treated as a scalar.
  818. c : array_like
  819. Array of coefficients ordered so that the coefficient of the term of
  820. multi-degree i,j,k is contained in ``c[i,j,k]``. If `c` has dimension
  821. greater than 3 the remaining indices enumerate multiple sets of
  822. coefficients.
  823. Returns
  824. -------
  825. values : ndarray, compatible object
  826. The values of the multidimensional polynomial on points formed with
  827. triples of corresponding values from `x`, `y`, and `z`.
  828. See Also
  829. --------
  830. hermeval, hermeval2d, hermegrid2d, hermegrid3d
  831. Notes
  832. -----
  833. .. versionadded:: 1.7.0
  834. """
  835. return pu._valnd(hermeval, c, x, y, z)
  836. def hermegrid3d(x, y, z, c):
  837. """
  838. Evaluate a 3-D HermiteE series on the Cartesian product of x, y, and z.
  839. This function returns the values:
  840. .. math:: p(a,b,c) = \\sum_{i,j,k} c_{i,j,k} * He_i(a) * He_j(b) * He_k(c)
  841. where the points `(a, b, c)` consist of all triples formed by taking
  842. `a` from `x`, `b` from `y`, and `c` from `z`. The resulting points form
  843. a grid with `x` in the first dimension, `y` in the second, and `z` in
  844. the third.
  845. The parameters `x`, `y`, and `z` are converted to arrays only if they
  846. are tuples or a lists, otherwise they are treated as a scalars. In
  847. either case, either `x`, `y`, and `z` or their elements must support
  848. multiplication and addition both with themselves and with the elements
  849. of `c`.
  850. If `c` has fewer than three dimensions, ones are implicitly appended to
  851. its shape to make it 3-D. The shape of the result will be c.shape[3:] +
  852. x.shape + y.shape + z.shape.
  853. Parameters
  854. ----------
  855. x, y, z : array_like, compatible objects
  856. The three dimensional series is evaluated at the points in the
  857. Cartesian product of `x`, `y`, and `z`. If `x`,`y`, or `z` is a
  858. list or tuple, it is first converted to an ndarray, otherwise it is
  859. left unchanged and, if it isn't an ndarray, it is treated as a
  860. scalar.
  861. c : array_like
  862. Array of coefficients ordered so that the coefficients for terms of
  863. degree i,j are contained in ``c[i,j]``. If `c` has dimension
  864. greater than two the remaining indices enumerate multiple sets of
  865. coefficients.
  866. Returns
  867. -------
  868. values : ndarray, compatible object
  869. The values of the two dimensional polynomial at points in the Cartesian
  870. product of `x` and `y`.
  871. See Also
  872. --------
  873. hermeval, hermeval2d, hermegrid2d, hermeval3d
  874. Notes
  875. -----
  876. .. versionadded:: 1.7.0
  877. """
  878. return pu._gridnd(hermeval, c, x, y, z)
  879. def hermevander(x, deg):
  880. """Pseudo-Vandermonde matrix of given degree.
  881. Returns the pseudo-Vandermonde matrix of degree `deg` and sample points
  882. `x`. The pseudo-Vandermonde matrix is defined by
  883. .. math:: V[..., i] = He_i(x),
  884. where `0 <= i <= deg`. The leading indices of `V` index the elements of
  885. `x` and the last index is the degree of the HermiteE polynomial.
  886. If `c` is a 1-D array of coefficients of length `n + 1` and `V` is the
  887. array ``V = hermevander(x, n)``, then ``np.dot(V, c)`` and
  888. ``hermeval(x, c)`` are the same up to roundoff. This equivalence is
  889. useful both for least squares fitting and for the evaluation of a large
  890. number of HermiteE series of the same degree and sample points.
  891. Parameters
  892. ----------
  893. x : array_like
  894. Array of points. The dtype is converted to float64 or complex128
  895. depending on whether any of the elements are complex. If `x` is
  896. scalar it is converted to a 1-D array.
  897. deg : int
  898. Degree of the resulting matrix.
  899. Returns
  900. -------
  901. vander : ndarray
  902. The pseudo-Vandermonde matrix. The shape of the returned matrix is
  903. ``x.shape + (deg + 1,)``, where The last index is the degree of the
  904. corresponding HermiteE polynomial. The dtype will be the same as
  905. the converted `x`.
  906. Examples
  907. --------
  908. >>> from numpy.polynomial.hermite_e import hermevander
  909. >>> x = np.array([-1, 0, 1])
  910. >>> hermevander(x, 3)
  911. array([[ 1., -1., 0., 2.],
  912. [ 1., 0., -1., -0.],
  913. [ 1., 1., 0., -2.]])
  914. """
  915. ideg = pu._deprecate_as_int(deg, "deg")
  916. if ideg < 0:
  917. raise ValueError("deg must be non-negative")
  918. x = np.array(x, copy=False, ndmin=1) + 0.0
  919. dims = (ideg + 1,) + x.shape
  920. dtyp = x.dtype
  921. v = np.empty(dims, dtype=dtyp)
  922. v[0] = x*0 + 1
  923. if ideg > 0:
  924. v[1] = x
  925. for i in range(2, ideg + 1):
  926. v[i] = (v[i-1]*x - v[i-2]*(i - 1))
  927. return np.moveaxis(v, 0, -1)
  928. def hermevander2d(x, y, deg):
  929. """Pseudo-Vandermonde matrix of given degrees.
  930. Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
  931. points `(x, y)`. The pseudo-Vandermonde matrix is defined by
  932. .. math:: V[..., (deg[1] + 1)*i + j] = He_i(x) * He_j(y),
  933. where `0 <= i <= deg[0]` and `0 <= j <= deg[1]`. The leading indices of
  934. `V` index the points `(x, y)` and the last index encodes the degrees of
  935. the HermiteE polynomials.
  936. If ``V = hermevander2d(x, y, [xdeg, ydeg])``, then the columns of `V`
  937. correspond to the elements of a 2-D coefficient array `c` of shape
  938. (xdeg + 1, ydeg + 1) in the order
  939. .. math:: c_{00}, c_{01}, c_{02} ... , c_{10}, c_{11}, c_{12} ...
  940. and ``np.dot(V, c.flat)`` and ``hermeval2d(x, y, c)`` will be the same
  941. up to roundoff. This equivalence is useful both for least squares
  942. fitting and for the evaluation of a large number of 2-D HermiteE
  943. series of the same degrees and sample points.
  944. Parameters
  945. ----------
  946. x, y : array_like
  947. Arrays of point coordinates, all of the same shape. The dtypes
  948. will be converted to either float64 or complex128 depending on
  949. whether any of the elements are complex. Scalars are converted to
  950. 1-D arrays.
  951. deg : list of ints
  952. List of maximum degrees of the form [x_deg, y_deg].
  953. Returns
  954. -------
  955. vander2d : ndarray
  956. The shape of the returned matrix is ``x.shape + (order,)``, where
  957. :math:`order = (deg[0]+1)*(deg[1]+1)`. The dtype will be the same
  958. as the converted `x` and `y`.
  959. See Also
  960. --------
  961. hermevander, hermevander3d, hermeval2d, hermeval3d
  962. Notes
  963. -----
  964. .. versionadded:: 1.7.0
  965. """
  966. return pu._vander_nd_flat((hermevander, hermevander), (x, y), deg)
  967. def hermevander3d(x, y, z, deg):
  968. """Pseudo-Vandermonde matrix of given degrees.
  969. Returns the pseudo-Vandermonde matrix of degrees `deg` and sample
  970. points `(x, y, z)`. If `l, m, n` are the given degrees in `x, y, z`,
  971. then Hehe pseudo-Vandermonde matrix is defined by
  972. .. math:: V[..., (m+1)(n+1)i + (n+1)j + k] = He_i(x)*He_j(y)*He_k(z),
  973. where `0 <= i <= l`, `0 <= j <= m`, and `0 <= j <= n`. The leading
  974. indices of `V` index the points `(x, y, z)` and the last index encodes
  975. the degrees of the HermiteE polynomials.
  976. If ``V = hermevander3d(x, y, z, [xdeg, ydeg, zdeg])``, then the columns
  977. of `V` correspond to the elements of a 3-D coefficient array `c` of
  978. shape (xdeg + 1, ydeg + 1, zdeg + 1) in the order
  979. .. math:: c_{000}, c_{001}, c_{002},... , c_{010}, c_{011}, c_{012},...
  980. and ``np.dot(V, c.flat)`` and ``hermeval3d(x, y, z, c)`` will be the
  981. same up to roundoff. This equivalence is useful both for least squares
  982. fitting and for the evaluation of a large number of 3-D HermiteE
  983. series of the same degrees and sample points.
  984. Parameters
  985. ----------
  986. x, y, z : array_like
  987. Arrays of point coordinates, all of the same shape. The dtypes will
  988. be converted to either float64 or complex128 depending on whether
  989. any of the elements are complex. Scalars are converted to 1-D
  990. arrays.
  991. deg : list of ints
  992. List of maximum degrees of the form [x_deg, y_deg, z_deg].
  993. Returns
  994. -------
  995. vander3d : ndarray
  996. The shape of the returned matrix is ``x.shape + (order,)``, where
  997. :math:`order = (deg[0]+1)*(deg[1]+1)*(deg[2]+1)`. The dtype will
  998. be the same as the converted `x`, `y`, and `z`.
  999. See Also
  1000. --------
  1001. hermevander, hermevander3d, hermeval2d, hermeval3d
  1002. Notes
  1003. -----
  1004. .. versionadded:: 1.7.0
  1005. """
  1006. return pu._vander_nd_flat((hermevander, hermevander, hermevander), (x, y, z), deg)
  1007. def hermefit(x, y, deg, rcond=None, full=False, w=None):
  1008. """
  1009. Least squares fit of Hermite series to data.
  1010. Return the coefficients of a HermiteE series of degree `deg` that is
  1011. the least squares fit to the data values `y` given at points `x`. If
  1012. `y` is 1-D the returned coefficients will also be 1-D. If `y` is 2-D
  1013. multiple fits are done, one for each column of `y`, and the resulting
  1014. coefficients are stored in the corresponding columns of a 2-D return.
  1015. The fitted polynomial(s) are in the form
  1016. .. math:: p(x) = c_0 + c_1 * He_1(x) + ... + c_n * He_n(x),
  1017. where `n` is `deg`.
  1018. Parameters
  1019. ----------
  1020. x : array_like, shape (M,)
  1021. x-coordinates of the M sample points ``(x[i], y[i])``.
  1022. y : array_like, shape (M,) or (M, K)
  1023. y-coordinates of the sample points. Several data sets of sample
  1024. points sharing the same x-coordinates can be fitted at once by
  1025. passing in a 2D-array that contains one dataset per column.
  1026. deg : int or 1-D array_like
  1027. Degree(s) of the fitting polynomials. If `deg` is a single integer
  1028. all terms up to and including the `deg`'th term are included in the
  1029. fit. For NumPy versions >= 1.11.0 a list of integers specifying the
  1030. degrees of the terms to include may be used instead.
  1031. rcond : float, optional
  1032. Relative condition number of the fit. Singular values smaller than
  1033. this relative to the largest singular value will be ignored. The
  1034. default value is len(x)*eps, where eps is the relative precision of
  1035. the float type, about 2e-16 in most cases.
  1036. full : bool, optional
  1037. Switch determining nature of return value. When it is False (the
  1038. default) just the coefficients are returned, when True diagnostic
  1039. information from the singular value decomposition is also returned.
  1040. w : array_like, shape (`M`,), optional
  1041. Weights. If not None, the contribution of each point
  1042. ``(x[i],y[i])`` to the fit is weighted by `w[i]`. Ideally the
  1043. weights are chosen so that the errors of the products ``w[i]*y[i]``
  1044. all have the same variance. The default value is None.
  1045. Returns
  1046. -------
  1047. coef : ndarray, shape (M,) or (M, K)
  1048. Hermite coefficients ordered from low to high. If `y` was 2-D,
  1049. the coefficients for the data in column k of `y` are in column
  1050. `k`.
  1051. [residuals, rank, singular_values, rcond] : list
  1052. These values are only returned if `full` = True
  1053. resid -- sum of squared residuals of the least squares fit
  1054. rank -- the numerical rank of the scaled Vandermonde matrix
  1055. sv -- singular values of the scaled Vandermonde matrix
  1056. rcond -- value of `rcond`.
  1057. For more details, see `linalg.lstsq`.
  1058. Warns
  1059. -----
  1060. RankWarning
  1061. The rank of the coefficient matrix in the least-squares fit is
  1062. deficient. The warning is only raised if `full` = False. The
  1063. warnings can be turned off by
  1064. >>> import warnings
  1065. >>> warnings.simplefilter('ignore', np.RankWarning)
  1066. See Also
  1067. --------
  1068. chebfit, legfit, polyfit, hermfit, polyfit
  1069. hermeval : Evaluates a Hermite series.
  1070. hermevander : pseudo Vandermonde matrix of Hermite series.
  1071. hermeweight : HermiteE weight function.
  1072. linalg.lstsq : Computes a least-squares fit from the matrix.
  1073. scipy.interpolate.UnivariateSpline : Computes spline fits.
  1074. Notes
  1075. -----
  1076. The solution is the coefficients of the HermiteE series `p` that
  1077. minimizes the sum of the weighted squared errors
  1078. .. math:: E = \\sum_j w_j^2 * |y_j - p(x_j)|^2,
  1079. where the :math:`w_j` are the weights. This problem is solved by
  1080. setting up the (typically) overdetermined matrix equation
  1081. .. math:: V(x) * c = w * y,
  1082. where `V` is the pseudo Vandermonde matrix of `x`, the elements of `c`
  1083. are the coefficients to be solved for, and the elements of `y` are the
  1084. observed values. This equation is then solved using the singular value
  1085. decomposition of `V`.
  1086. If some of the singular values of `V` are so small that they are
  1087. neglected, then a `RankWarning` will be issued. This means that the
  1088. coefficient values may be poorly determined. Using a lower order fit
  1089. will usually get rid of the warning. The `rcond` parameter can also be
  1090. set to a value smaller than its default, but the resulting fit may be
  1091. spurious and have large contributions from roundoff error.
  1092. Fits using HermiteE series are probably most useful when the data can
  1093. be approximated by ``sqrt(w(x)) * p(x)``, where `w(x)` is the HermiteE
  1094. weight. In that case the weight ``sqrt(w(x[i]))`` should be used
  1095. together with data values ``y[i]/sqrt(w(x[i]))``. The weight function is
  1096. available as `hermeweight`.
  1097. References
  1098. ----------
  1099. .. [1] Wikipedia, "Curve fitting",
  1100. https://en.wikipedia.org/wiki/Curve_fitting
  1101. Examples
  1102. --------
  1103. >>> from numpy.polynomial.hermite_e import hermefit, hermeval
  1104. >>> x = np.linspace(-10, 10)
  1105. >>> np.random.seed(123)
  1106. >>> err = np.random.randn(len(x))/10
  1107. >>> y = hermeval(x, [1, 2, 3]) + err
  1108. >>> hermefit(x, y, 2)
  1109. array([ 1.01690445, 1.99951418, 2.99948696]) # may vary
  1110. """
  1111. return pu._fit(hermevander, x, y, deg, rcond, full, w)
  1112. def hermecompanion(c):
  1113. """
  1114. Return the scaled companion matrix of c.
  1115. The basis polynomials are scaled so that the companion matrix is
  1116. symmetric when `c` is an HermiteE basis polynomial. This provides
  1117. better eigenvalue estimates than the unscaled case and for basis
  1118. polynomials the eigenvalues are guaranteed to be real if
  1119. `numpy.linalg.eigvalsh` is used to obtain them.
  1120. Parameters
  1121. ----------
  1122. c : array_like
  1123. 1-D array of HermiteE series coefficients ordered from low to high
  1124. degree.
  1125. Returns
  1126. -------
  1127. mat : ndarray
  1128. Scaled companion matrix of dimensions (deg, deg).
  1129. Notes
  1130. -----
  1131. .. versionadded:: 1.7.0
  1132. """
  1133. # c is a trimmed copy
  1134. [c] = pu.as_series([c])
  1135. if len(c) < 2:
  1136. raise ValueError('Series must have maximum degree of at least 1.')
  1137. if len(c) == 2:
  1138. return np.array([[-c[0]/c[1]]])
  1139. n = len(c) - 1
  1140. mat = np.zeros((n, n), dtype=c.dtype)
  1141. scl = np.hstack((1., 1./np.sqrt(np.arange(n - 1, 0, -1))))
  1142. scl = np.multiply.accumulate(scl)[::-1]
  1143. top = mat.reshape(-1)[1::n+1]
  1144. bot = mat.reshape(-1)[n::n+1]
  1145. top[...] = np.sqrt(np.arange(1, n))
  1146. bot[...] = top
  1147. mat[:, -1] -= scl*c[:-1]/c[-1]
  1148. return mat
  1149. def hermeroots(c):
  1150. """
  1151. Compute the roots of a HermiteE series.
  1152. Return the roots (a.k.a. "zeros") of the polynomial
  1153. .. math:: p(x) = \\sum_i c[i] * He_i(x).
  1154. Parameters
  1155. ----------
  1156. c : 1-D array_like
  1157. 1-D array of coefficients.
  1158. Returns
  1159. -------
  1160. out : ndarray
  1161. Array of the roots of the series. If all the roots are real,
  1162. then `out` is also real, otherwise it is complex.
  1163. See Also
  1164. --------
  1165. polyroots, legroots, lagroots, hermroots, chebroots
  1166. Notes
  1167. -----
  1168. The root estimates are obtained as the eigenvalues of the companion
  1169. matrix, Roots far from the origin of the complex plane may have large
  1170. errors due to the numerical instability of the series for such
  1171. values. Roots with multiplicity greater than 1 will also show larger
  1172. errors as the value of the series near such points is relatively
  1173. insensitive to errors in the roots. Isolated roots near the origin can
  1174. be improved by a few iterations of Newton's method.
  1175. The HermiteE series basis polynomials aren't powers of `x` so the
  1176. results of this function may seem unintuitive.
  1177. Examples
  1178. --------
  1179. >>> from numpy.polynomial.hermite_e import hermeroots, hermefromroots
  1180. >>> coef = hermefromroots([-1, 0, 1])
  1181. >>> coef
  1182. array([0., 2., 0., 1.])
  1183. >>> hermeroots(coef)
  1184. array([-1., 0., 1.]) # may vary
  1185. """
  1186. # c is a trimmed copy
  1187. [c] = pu.as_series([c])
  1188. if len(c) <= 1:
  1189. return np.array([], dtype=c.dtype)
  1190. if len(c) == 2:
  1191. return np.array([-c[0]/c[1]])
  1192. # rotated companion matrix reduces error
  1193. m = hermecompanion(c)[::-1,::-1]
  1194. r = la.eigvals(m)
  1195. r.sort()
  1196. return r
  1197. def _normed_hermite_e_n(x, n):
  1198. """
  1199. Evaluate a normalized HermiteE polynomial.
  1200. Compute the value of the normalized HermiteE polynomial of degree ``n``
  1201. at the points ``x``.
  1202. Parameters
  1203. ----------
  1204. x : ndarray of double.
  1205. Points at which to evaluate the function
  1206. n : int
  1207. Degree of the normalized HermiteE function to be evaluated.
  1208. Returns
  1209. -------
  1210. values : ndarray
  1211. The shape of the return value is described above.
  1212. Notes
  1213. -----
  1214. .. versionadded:: 1.10.0
  1215. This function is needed for finding the Gauss points and integration
  1216. weights for high degrees. The values of the standard HermiteE functions
  1217. overflow when n >= 207.
  1218. """
  1219. if n == 0:
  1220. return np.full(x.shape, 1/np.sqrt(np.sqrt(2*np.pi)))
  1221. c0 = 0.
  1222. c1 = 1./np.sqrt(np.sqrt(2*np.pi))
  1223. nd = float(n)
  1224. for i in range(n - 1):
  1225. tmp = c0
  1226. c0 = -c1*np.sqrt((nd - 1.)/nd)
  1227. c1 = tmp + c1*x*np.sqrt(1./nd)
  1228. nd = nd - 1.0
  1229. return c0 + c1*x
  1230. def hermegauss(deg):
  1231. """
  1232. Gauss-HermiteE quadrature.
  1233. Computes the sample points and weights for Gauss-HermiteE quadrature.
  1234. These sample points and weights will correctly integrate polynomials of
  1235. degree :math:`2*deg - 1` or less over the interval :math:`[-\\inf, \\inf]`
  1236. with the weight function :math:`f(x) = \\exp(-x^2/2)`.
  1237. Parameters
  1238. ----------
  1239. deg : int
  1240. Number of sample points and weights. It must be >= 1.
  1241. Returns
  1242. -------
  1243. x : ndarray
  1244. 1-D ndarray containing the sample points.
  1245. y : ndarray
  1246. 1-D ndarray containing the weights.
  1247. Notes
  1248. -----
  1249. .. versionadded:: 1.7.0
  1250. The results have only been tested up to degree 100, higher degrees may
  1251. be problematic. The weights are determined by using the fact that
  1252. .. math:: w_k = c / (He'_n(x_k) * He_{n-1}(x_k))
  1253. where :math:`c` is a constant independent of :math:`k` and :math:`x_k`
  1254. is the k'th root of :math:`He_n`, and then scaling the results to get
  1255. the right value when integrating 1.
  1256. """
  1257. ideg = pu._deprecate_as_int(deg, "deg")
  1258. if ideg <= 0:
  1259. raise ValueError("deg must be a positive integer")
  1260. # first approximation of roots. We use the fact that the companion
  1261. # matrix is symmetric in this case in order to obtain better zeros.
  1262. c = np.array([0]*deg + [1])
  1263. m = hermecompanion(c)
  1264. x = la.eigvalsh(m)
  1265. # improve roots by one application of Newton
  1266. dy = _normed_hermite_e_n(x, ideg)
  1267. df = _normed_hermite_e_n(x, ideg - 1) * np.sqrt(ideg)
  1268. x -= dy/df
  1269. # compute the weights. We scale the factor to avoid possible numerical
  1270. # overflow.
  1271. fm = _normed_hermite_e_n(x, ideg - 1)
  1272. fm /= np.abs(fm).max()
  1273. w = 1/(fm * fm)
  1274. # for Hermite_e we can also symmetrize
  1275. w = (w + w[::-1])/2
  1276. x = (x - x[::-1])/2
  1277. # scale w to get the right value
  1278. w *= np.sqrt(2*np.pi) / w.sum()
  1279. return x, w
  1280. def hermeweight(x):
  1281. """Weight function of the Hermite_e polynomials.
  1282. The weight function is :math:`\\exp(-x^2/2)` and the interval of
  1283. integration is :math:`[-\\inf, \\inf]`. the HermiteE polynomials are
  1284. orthogonal, but not normalized, with respect to this weight function.
  1285. Parameters
  1286. ----------
  1287. x : array_like
  1288. Values at which the weight function will be computed.
  1289. Returns
  1290. -------
  1291. w : ndarray
  1292. The weight function at `x`.
  1293. Notes
  1294. -----
  1295. .. versionadded:: 1.7.0
  1296. """
  1297. w = np.exp(-.5*x**2)
  1298. return w
  1299. #
  1300. # HermiteE series class
  1301. #
  1302. class HermiteE(ABCPolyBase):
  1303. """An HermiteE series class.
  1304. The HermiteE class provides the standard Python numerical methods
  1305. '+', '-', '*', '//', '%', 'divmod', '**', and '()' as well as the
  1306. attributes and methods listed in the `ABCPolyBase` documentation.
  1307. Parameters
  1308. ----------
  1309. coef : array_like
  1310. HermiteE coefficients in order of increasing degree, i.e,
  1311. ``(1, 2, 3)`` gives ``1*He_0(x) + 2*He_1(X) + 3*He_2(x)``.
  1312. domain : (2,) array_like, optional
  1313. Domain to use. The interval ``[domain[0], domain[1]]`` is mapped
  1314. to the interval ``[window[0], window[1]]`` by shifting and scaling.
  1315. The default value is [-1, 1].
  1316. window : (2,) array_like, optional
  1317. Window, see `domain` for its use. The default value is [-1, 1].
  1318. .. versionadded:: 1.6.0
  1319. """
  1320. # Virtual Functions
  1321. _add = staticmethod(hermeadd)
  1322. _sub = staticmethod(hermesub)
  1323. _mul = staticmethod(hermemul)
  1324. _div = staticmethod(hermediv)
  1325. _pow = staticmethod(hermepow)
  1326. _val = staticmethod(hermeval)
  1327. _int = staticmethod(hermeint)
  1328. _der = staticmethod(hermeder)
  1329. _fit = staticmethod(hermefit)
  1330. _line = staticmethod(hermeline)
  1331. _roots = staticmethod(hermeroots)
  1332. _fromroots = staticmethod(hermefromroots)
  1333. # Virtual properties
  1334. nickname = 'herme'
  1335. domain = np.array(hermedomain)
  1336. window = np.array(hermedomain)
  1337. basis_name = 'He'