gitee.com/quant1x/num@v0.3.2/polynomial.go (about) 1 package num 2 3 import ( 4 "gitee.com/quant1x/num/x64" 5 "math" 6 ) 7 8 // PolyFit 9 // 10 // Least squares polynomial fit. 11 // 12 // .. note:: 13 // This forms part of the old polynomial API. Since version 1.4, the 14 // new polynomial API defined in `numpy.polynomial` is preferred. 15 // A summary of the differences can be found in the 16 // :doc:`transition guide </reference/routines.polynomials>`. 17 // 18 // Fit a polynomial ``p(x) = p[0] * x**deg + ... + p[deg]`` of degree `deg` 19 // to points `(x, y)`. Returns a vector of coefficients `p` that minimises 20 // the squared error in the order `deg`, `deg-1`, ... `0`. 21 // 22 // The `Polynomial.fit <numpy.polynomial.polynomial.Polynomial.fit>` class 23 // method is recommended for new code as it is more stable numerically. See 24 // the documentation of the method for more information. 25 // 26 // Parameters 27 // ---------- 28 // x : array_like, shape (M,) 29 // x-coordinates of the M sample points ``(x[i], y[i])``. 30 // y : array_like, shape (M,) or (M, K) 31 // y-coordinates of the sample points. Several data sets of sample 32 // points sharing the same x-coordinates can be fitted at once by 33 // passing in a 2D-array that contains one dataset per column. 34 // deg : int 35 // Degree of the fitting polynomial 36 // 37 // Returns 38 // ------- 39 // p : ndarray, shape (deg + 1,) or (deg + 1, K) 40 // Polynomial coefficients, highest power first. If `y` was 2-D, the 41 // coefficients for `k`-th data set are in ``p[:,k]``. 42 // 43 // residuals, rank, singular_values, rcond 44 // These values are only returned if ``full == True`` 45 // 46 // - residuals -- sum of squared residuals of the least squares fit 47 // - rank -- the effective rank of the scaled Vandermonde 48 // coefficient matrix 49 // - singular_values -- singular values of the scaled Vandermonde 50 // coefficient matrix 51 // - rcond -- value of `rcond`. 52 // 53 // For more details, see `numpy.linalg.lstsq`. 54 // 55 // Warns 56 // ----- 57 // RankWarning 58 // The rank of the coefficient matrix in the least-squares fit is 59 // deficient. The warning is only raised if ``full == False``. 60 // 61 // The warnings can be turned off by 62 // 63 // >>> import warnings 64 // >>> warnings.simplefilter('ignore', np.RankWarning) 65 // 66 // See Also 67 // -------- 68 // polyval : Compute polynomial values. 69 // linalg.lstsq : Computes a least-squares fit. 70 // scipy.interpolate.UnivariateSpline : Computes spline fits. 71 // 72 // Notes 73 // ----- 74 // The solution minimizes the squared error 75 // 76 // .. math:: 77 // E = \\sum_{j=0}^k |p(x_j) - y_j|^2 78 // 79 // in the equations:: 80 // 81 // x[0]**n * p[0] + ... + x[0] * p[n-1] + p[n] = y[0] 82 // x[1]**n * p[0] + ... + x[1] * p[n-1] + p[n] = y[1] 83 // ... 84 // x[k]**n * p[0] + ... + x[k] * p[n-1] + p[n] = y[k] 85 // 86 // The coefficient matrix of the coefficients `p` is a Vandermonde matrix. 87 // 88 // `polyfit` issues a `RankWarning` when the least-squares fit is badly 89 // conditioned. This implies that the best fit is not well-defined due 90 // to numerical error. The results may be improved by lowering the polynomial 91 // degree or by replacing `x` by `x` - `x`.mean(). The `rcond` parameter 92 // can also be set to a value smaller than its default, but the resulting 93 // fit may be spurious: including contributions from the small singular 94 // values can add numerical noise to the result. 95 // 96 // Note that fitting polynomial coefficients is inherently badly conditioned 97 // when the degree of the polynomial is large or the interval of sample points 98 // is badly centered. The quality of the fit should always be checked in these 99 // cases. When polynomial fits are not satisfactory, splines may be a good 100 // alternative. 101 // 102 // References 103 // ---------- 104 // .. [1] Wikipedia, "Curve fitting", 105 // https://en.wikipedia.org/wiki/Curve_fitting 106 // .. [2] Wikipedia, "Polynomial interpolation", 107 // https://en.wikipedia.org/wiki/Polynomial_interpolation 108 // .. [3] numpy.polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False) 109 // https://numpy.org/doc/stable/reference/generated/numpy.polyfit.html 110 // 111 // Examples 112 // -------- 113 // >>> import warnings 114 // >>> x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0]) 115 // >>> y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0]) 116 // >>> z = np.polyfit(x, y, 3) 117 // >>> z 118 // array([ 0.08703704, -0.81349206, 1.69312169, -0.03968254]) # may vary 119 func PolyFit(x, y []DType, deg int, args ...any) []DType { 120 // 默认从右向左 121 var __increasing = false 122 if len(args) > 0 { 123 // 第一个参数为是否copy 124 if arg0, ok := args[0].(bool); ok { 125 __increasing = arg0 126 } 127 } 128 // Initialize matrix to store powers of x 129 var X = make([][]float64, len(x)) 130 for i := range X { 131 X[i] = make([]float64, deg+1) 132 } 133 134 // Fill matrix with powers of x 135 for i := 0; i < len(x); i++ { 136 for j := 0; j <= deg; j++ { 137 k := j 138 if !__increasing { 139 k = deg - k 140 } 141 X[i][j] = math.Pow(x[i], float64(k)) 142 } 143 } 144 145 // Calculate transpose of X 146 var XT = __transpose(X) 147 148 // Multiply XT with X 149 var XTX = __matmul(XT, X) 150 151 // Invert XTX 152 var XTXinv = __inv(XTX) 153 154 // Multiply XTXinv with XT 155 var XTXinvXT = __matmul(XTXinv, XT) 156 157 // Multiply result with y 158 var coef = __matvec(XTXinvXT, y) 159 160 return coef 161 } 162 163 // PolyVal 164 // 165 // Evaluate a polynomial at specific values. 166 // 167 // .. note:: 168 // This forms part of the old polynomial API. Since version 1.4, the 169 // new polynomial API defined in `numpy.polynomial` is preferred. 170 // A summary of the differences can be found in the 171 // :doc:`transition guide </reference/routines.polynomials>`. 172 // 173 // If `p` is of length N, this function returns the value: 174 // 175 // ``p[0]*x**(N-1) + p[1]*x**(N-2) + ... + p[N-2]*x + p[N-1]`` 176 // 177 // If `x` is a sequence, then ``p(x)`` is returned for each element of ``x``. 178 // If `x` is another polynomial then the composite polynomial ``p(x(t))`` 179 // is returned. 180 // 181 // Parameters 182 // 183 // ---------- 184 // 185 // p : array_like or poly1d object 186 // 1D array of polynomial coefficients (including coefficients equal 187 // to zero) from highest degree to the constant term, or an 188 // instance of poly1d. 189 // x : array_like or poly1d object 190 // A number, an array of numbers, or an instance of poly1d, at 191 // which to evaluate `p`. 192 // 193 // Returns 194 // 195 // ------- 196 // 197 // values : ndarray or poly1d 198 // If `x` is a poly1d instance, the result is the composition of the two 199 // polynomials, i.e., `x` is "substituted" in `p` and the simplified 200 // result is returned. In addition, the type of `x` - array_like or 201 // poly1d - governs the type of the output: `x` array_like => `values` 202 // array_like, `x` a poly1d object => `values` is also. 203 // 204 // See Also 205 // -------- 206 // poly1d: A polynomial class. 207 // 208 // Notes 209 // ----- 210 // Horner's scheme [1]_ is used to evaluate the polynomial. Even so, 211 // for polynomials of high degree the values may be inaccurate due to 212 // rounding errors. Use carefully. 213 // 214 // If `x` is a subtype of `ndarray` the return value will be of the same type. 215 // 216 // References 217 // ---------- 218 // .. [1] I. N. Bronshtein, K. A. Semendyayev, and K. A. Hirsch (Eng. 219 // trans. Ed.), *Handbook of Mathematics*, New York, Van Nostrand 220 // Reinhold Co., 1985, pg. 720. 221 // 222 // Examples 223 // -------- 224 // >>> np.polyval([3,0,1], 5) # 3 * 5**2 + 0 * 5**1 + 1 225 // 76 226 // >>> np.polyval([3,0,1], np.poly1d(5)) 227 // poly1d([76]) 228 // >>> np.polyval(np.poly1d([3,0,1]), 5) 229 // 76 230 // >>> np.polyval(np.poly1d([3,0,1]), np.poly1d(5)) 231 // poly1d([76]) 232 func PolyVal(p, x []DType) []DType { 233 //p = NX.asarray(p) 234 //if isinstance(x, poly1d): 235 // y = 0 236 //else: 237 // x = NX.asanyarray(x) 238 // y = NX.zeros_like(x) 239 //for pv in p: 240 //y = y * x + pv 241 y := Repeat(DType(0), len(x)) 242 for _, v := range p { 243 x64.Mul_Inplace(y, x) 244 x64.AddNumber_Inplace(y, v) 245 } 246 return y 247 } 248 249 // Function to calculate matrix transpose 250 func __transpose(a [][]float64) [][]float64 { 251 var transposed = make([][]float64, len(a[0])) 252 for i := range transposed { 253 transposed[i] = make([]float64, len(a)) 254 } 255 for i := 0; i < len(a); i++ { 256 for j := 0; j < len(a[0]); j++ { 257 transposed[j][i] = a[i][j] 258 } 259 } 260 return transposed 261 } 262 263 // Function to multiply two matrices 264 func __matmul(a, b [][]float64) [][]float64 { 265 var result = make([][]float64, len(a)) 266 for i := range result { 267 result[i] = make([]float64, len(b[0])) 268 } 269 for i := 0; i < len(a); i++ { 270 for j := 0; j < len(b[0]); j++ { 271 for k := 0; k < len(b); k++ { 272 result[i][j] += a[i][k] * b[k][j] 273 } 274 } 275 } 276 return result 277 } 278 279 // Function to multiply a matrix with a vector 280 func __matvec(a [][]float64, b []float64) []float64 { 281 var result = make([]float64, len(a)) 282 for i := 0; i < len(a); i++ { 283 for j := 0; j < len(b); j++ { 284 result[i] += a[i][j] * b[j] 285 } 286 } 287 return result 288 } 289 290 func __inv(a [][]float64) [][]float64 { 291 var n = len(a) 292 293 // Create augmented matrix 294 var augmented = make([][]float64, n) 295 for i := range augmented { 296 augmented[i] = make([]float64, 2*n) 297 for j := 0; j < n; j++ { 298 augmented[i][j] = a[i][j] 299 } 300 } 301 for i := 0; i < n; i++ { 302 augmented[i][i+n] = 1 303 } 304 305 // Perform Gaussian elimination 306 for i := 0; i < n; i++ { 307 var pivot = augmented[i][i] 308 for j := i + 1; j < n; j++ { 309 var factor = augmented[j][i] / pivot 310 for k := i; k < 2*n; k++ { 311 augmented[j][k] -= factor * augmented[i][k] 312 } 313 } 314 } 315 316 // Perform back-substitution 317 for i := n - 1; i >= 0; i-- { 318 var pivot = augmented[i][i] 319 for j := i - 1; j >= 0; j-- { 320 var factor = augmented[j][i] / pivot 321 for k := i; k < 2*n; k++ { 322 augmented[j][k] -= factor * augmented[i][k] 323 } 324 } 325 } 326 327 // Normalize rows 328 for i := 0; i < n; i++ { 329 var pivot = augmented[i][i] 330 for j := 0; j < 2*n; j++ { 331 augmented[i][j] /= pivot 332 } 333 } 334 335 // Extract inverse from augmented matrix 336 var inverse = make([][]float64, n) 337 for i := range inverse { 338 inverse[i] = make([]float64, n) 339 } 340 for i := 0; i < n; i++ { 341 for j := 0; j < n; j++ { 342 inverse[i][j] = augmented[i][j+n] 343 } 344 } 345 346 return inverse 347 }