关于非线性曲线拟合的问题

rubyish

a hard way ~
f(x) = (17.800935828877)*x**2 + (7.19552718360075)*x

1. #!/usr/bin/perl
   
2. 
   
3. # queke:
   
4. # gnuplot> f(x) = (17.800935828877)*x**2 + (7.19552718360075)*x
   
5. # gnuplot> plot 'data' with points ps 3 pt 7, f(x) lw 5
   
6. 
   
7. use 5.016;
   
8. my ( $M, $Y );
   
9. 
   
10. while (<DATA>) {
    
11.     my ( $x, $y ) = split;
    
12.     push @$Y, [$y];
    
13.     push @$M, [ $x, $x * $x ];
    
14. }
    
15. 
    
16. sub gimme {
    
17.     my ( $M, $Y ) = @_;
    
18.     my $TM = T($M);
    
19.     my @A = map $_->[0], @{ M( I( M( $TM, $M ) ), M( $TM, $Y ) ) };
    
20.     "f(x) = ($A[1])*x**2 + ($A[0])*x";
    
21. }
    
22. 
    
23. my $fx = gimme $M, $Y;
    
24. say $fx;
    
25. 
    
26. sub T {
    
27.     my $m = shift;
    
28.     my @M = map {
    
29.         my $i = $_;
    
30.         [ map { $m->[$_][$i] } 0 .. $#$m ]
    
31.     } 0 .. $#{ $m->[0] };
    
32.     [@M];
    
33. }
    
34. 
    
35. sub M {
    
36.     my ( $A, $B ) = @_;
    
37.     my @M = map {
    
38.         my $i = $_;
    
39.         my @S = map {
    
40.             my ( $j, $S ) = $_;
    
41.             $S += $A->[$i][$_] * $B->[$_][$j] for 0 .. $#{$B};
    
42.             $S;
    
43.         } 0 .. $#{ $B->[0] };
    
44.         [@S]
    
45.     } 0 .. $#{$A};
    
46.     [@M];
    
47. }
    
48. 
    
49. sub I {
    
50.     my @m = @{ +shift };
    
51.     for my $i ( 0 .. $#m ) {
    
52.         my $it = 1.0 / $m[$i][$i];
    
53.         $m[$i][$i] = $it;
    
54.         for my $j ( 0 .. $#m ) {
    
55.             next if $j == $i;
    
56.             $m[$i][$j] *= -$it;
    
57.             $m[$j][$i] *= $it;
    
58.         }
    
59.         for my $j ( 0 .. $#m ) {
    
60.             next if $j == $i;
    
61.             for my $k ( 0 .. $#m ) {
    
62.                 next if $k == $i;
    
63.                 $m[$j][$k] += $m[$j][$i] * $m[$i][$k] / $it;
    
64.             }
    
65.         }
    
66.     }
    
67.     [@m];
    
68. }
    
69. 
    
70. 
    
71. 
    
72. __DATA__
    
73. 0       0
    
74. 0.05    0.5487
    
75. 0.1     1.0967
    
76. 0.15    1.6447
    
77. 0.2     2.1937
    
78. 0.25    2.7417
    
79. 0.3     3.5587
    
80. 0.35    4.6387
    
81. 0.4     5.7407
    
82. 0.45    6.9717
    

复制代码

rubyish

回复 11# minirain


    my $a1 = shift;
    my @m  = @{$a1};