数式微分怖いです

ひゃー。

> eval(parse(text=paste("D(expression(",gauusian(x),"),\"h\")")))
-(exp(-1/2 * ((Xi - 25)/h)^2) * (-1/2 * (2 * ((Xi - 25)/h^2 * 
    ((Xi - 25)/h)))) + (exp(-1/2 * ((Xi - 24)/h)^2) * (-1/2 * 
    (2 * ((Xi - 24)/h^2 * ((Xi - 24)/h)))) + (exp(-1/2 * ((Xi - 
    24)/h)^2) * (-1/2 * (2 * ((Xi - 24)/h^2 * ((Xi - 24)/h)))) + 
    (exp(-1/2 * ((Xi - 24)/h)^2) * (-1/2 * (2 * ((Xi - 24)/h^2 * 
        ((Xi - 24)/h)))) + (exp(-1/2 * ((Xi - 24)/h)^2) * (-1/2 * 
        (2 * ((Xi - 24)/h^2 * ((Xi - 24)/h)))) + (exp(-1/2 * 
        ((Xi - 23)/h)^2) * (-1/2 * (2 * ((Xi - 23)/h^2 * ((Xi - 
        23)/h)))) + (exp(-1/2 * ((Xi - 22)/h)^2) * (-1/2 * (2 * 
        ((Xi - 22)/h^2 * ((Xi - 22)/h)))) + (exp(-1/2 * ((Xi - 
        20)/h)^2) * (-1/2 * (2 * ((Xi - 20)/h^2 * ((Xi - 20)/h)))) + 
        (exp(-1/2 * ((Xi - 20)/h)^2) * (-1/2 * (2 * ((Xi - 20)/h^2 * 
            ((Xi - 20)/h)))) + (exp(-1/2 * ((Xi - 20)/h)^2) * 
            (-1/2 * (2 * ((Xi - 20)/h^2 * ((Xi - 20)/h)))) + 
            (exp(-1/2 * ((Xi - 20)/h)^2) * (-1/2 * (2 * ((Xi - 
                20)/h^2 * ((Xi - 20)/h)))) + (exp(-1/2 * ((Xi - 
                20)/h)^2) * (-1/2 * (2 * ((Xi - 20)/h^2 * ((Xi - 
                20)/h)))) + (exp(-1/2 * ((Xi - 19)/h)^2) * (-1/2 * 
                (2 * ((Xi - 19)/h^2 * ((Xi - 19)/h)))) + (exp(-1/2 * 
                ((Xi - 19)/h)^2) * (-1/2 * (2 * ((Xi - 19)/h^2 * 
                ((Xi - 19)/h)))) + (exp(-1/2 * ((Xi - 19)/h)^2) * 
                (-1/2 * (2 * ((Xi - 19)/h^2 * ((Xi - 19)/h)))) + 
                (exp(-1/2 * ((Xi - 18)/h)^2) * (-1/2 * (2 * ((Xi - 
                  18)/h^2 * ((Xi - 18)/h)))) + (exp(-1/2 * ((Xi - 
                  18)/h)^2) * (-1/2 * (2 * ((Xi - 18)/h^2 * ((Xi - 
                  18)/h)))) + (exp(-1/2 * ((Xi - 18)/h)^2) * 
                  (-1/2 * (2 * ((Xi - 18)/h^2 * ((Xi - 18)/h)))) + 
                  (exp(-1/2 * ((Xi - 18)/h)^2) * (-1/2 * (2 * 
                    ((Xi - 18)/h^2 * ((Xi - 18)/h)))) + (exp(-1/2 * 
                    ((Xi - 17)/h)^2) * (-1/2 * (2 * ((Xi - 17)/h^2 * 
                    ((Xi - 17)/h)))) + (exp(-1/2 * ((Xi - 17)/h)^2) * 
                    (-1/2 * (2 * ((Xi - 17)/h^2 * ((Xi - 17)/h)))) + 
                    (exp(-1/2 * ((Xi - 17)/h)^2) * (-1/2 * (2 * 
                      ((Xi - 17)/h^2 * ((Xi - 17)/h)))) + (exp(-1/2 * 
                      ((Xi - 16)/h)^2) * (-1/2 * (2 * ((Xi - 
                      16)/h^2 * ((Xi - 16)/h)))) + (exp(-1/2 * 
                      ((Xi - 16)/h)^2) * (-1/2 * (2 * ((Xi - 
                      16)/h^2 * ((Xi - 16)/h)))) + (exp(-1/2 * 
                      ((Xi - 15)/h)^2) * (-1/2 * (2 * ((Xi - 
                      15)/h^2 * ((Xi - 15)/h)))) + (exp(-1/2 * 
                      ((Xi - 15)/h)^2) * (-1/2 * (2 * ((Xi - 
                      15)/h^2 * ((Xi - 15)/h)))) + (exp(-1/2 * 
                      ((Xi - 15)/h)^2) * (-1/2 * (2 * ((Xi - 
                      15)/h^2 * ((Xi - 15)/h)))) + (exp(-1/2 * 
                      ((Xi - 14)/h)^2) * (-1/2 * (2 * ((Xi - 
                      14)/h^2 * ((Xi - 14)/h)))) + (exp(-1/2 * 
                      ((Xi - 14)/h)^2) * (-1/2 * (2 * ((Xi - 
                      14)/h^2 * ((Xi - 14)/h)))) + (exp(-1/2 * 
                      ((Xi - 14)/h)^2) * (-1/2 * (2 * ((Xi - 
                      14)/h^2 * ((Xi - 14)/h)))) + (exp(-1/2 * 
                      ((Xi - 14)/h)^2) * (-1/2 * (2 * ((Xi - 
                      14)/h^2 * ((Xi - 14)/h)))) + (exp(-1/2 * 
                      ((Xi - 13)/h)^2) * (-1/2 * (2 * ((Xi - 
                      13)/h^2 * ((Xi - 13)/h)))) + (exp(-1/2 * 
                      ((Xi - 13)/h)^2) * (-1/2 * (2 * ((Xi - 
                      13)/h^2 * ((Xi - 13)/h)))) + (exp(-1/2 * 
                      ((Xi - 13)/h)^2) * (-1/2 * (2 * ((Xi - 
                      13)/h^2 * ((Xi - 13)/h)))) + (exp(-1/2 * 
                      ((Xi - 13)/h)^2) * (-1/2 * (2 * ((Xi - 
                      13)/h^2 * ((Xi - 13)/h)))) + (exp(-1/2 * 
                      ((Xi - 12)/h)^2) * (-1/2 * (2 * ((Xi - 
                      12)/h^2 * ((Xi - 12)/h)))) + (exp(-1/2 * 
                      ((Xi - 12)/h)^2) * (-1/2 * (2 * ((Xi - 
                      12)/h^2 * ((Xi - 12)/h)))) + (exp(-1/2 * 
                      ((Xi - 12)/h)^2) * (-1/2 * (2 * ((Xi - 
                      12)/h^2 * ((Xi - 12)/h)))) + (exp(-1/2 * 
                      ((Xi - 12)/h)^2) * (-1/2 * (2 * ((Xi - 
                      12)/h^2 * ((Xi - 12)/h)))) + (exp(-1/2 * 
                      ((Xi - 11)/h)^2) * (-1/2 * (2 * ((Xi - 
                      11)/h^2 * ((Xi - 11)/h)))) + (exp(-1/2 * 
                      ((Xi - 11)/h)^2) * (-1/2 * (2 * ((Xi - 
                      11)/h^2 * ((Xi - 11)/h)))) + (exp(-1/2 * 
                      ((Xi - 10)/h)^2) * (-1/2 * (2 * ((Xi - 
                      10)/h^2 * ((Xi - 10)/h)))) + (exp(-1/2 * 
                      ((Xi - 10)/h)^2) * (-1/2 * (2 * ((Xi - 
                      10)/h^2 * ((Xi - 10)/h)))) + (exp(-1/2 * 
                      ((Xi - 10)/h)^2) * (-1/2 * (2 * ((Xi - 
                      10)/h^2 * ((Xi - 10)/h)))) + (exp(-1/2 * 
                      ((Xi - 9)/h)^2) * (-1/2 * (2 * ((Xi - 9)/h^2 * 
                      ((Xi - 9)/h)))) + (exp(-1/2 * ((Xi - 8)/h)^2) * 
                      (-1/2 * (2 * ((Xi - 8)/h^2 * ((Xi - 8)/h)))) + 
                      (exp(-1/2 * ((Xi - 7)/h)^2) * (-1/2 * (2 * 
                        ((Xi - 7)/h^2 * ((Xi - 7)/h)))) + (exp(-1/2 * 
                        ((Xi - 7)/h)^2) * (-1/2 * (2 * ((Xi - 
                        7)/h^2 * ((Xi - 7)/h)))) + (exp(-1/2 * 
                        ((Xi - 4)/h)^2) * (-1/2 * (2 * ((Xi - 
                        4)/h^2 * ((Xi - 4)/h)))) + exp(-1/2 * 
                        ((Xi - 4)/h)^2) * (-1/2 * (2 * ((Xi - 
                        4)/h^2 * ((Xi - 4)/h)))))))))))))))))))))))))))))))))))))))))))))))))))))

この辺を使いつつ、以下の最小二乗クロスバリデーション評価関数を最小化します。

> cv(x,y)
[1] "( 2 - exp(-1/2 * ((Xi - 4)/h)^2) * 2  - Yi / exp(-1/2 * ((Xi - 4)/h)^2) - 1 + 10 - exp(-1/2 * ((Xi - 4)/h)^2) * 10  - Yi / exp(-1/2 * ((Xi - 4)/h)^2) - 1 + 4 - exp(-1/2 * ((Xi - 7)/h)^2) * 4  - Yi / exp(-1/2 * ((Xi - 7)/h)^2) - 1 + 22 - exp(-1/2 * ((Xi - 7)/h)^2) * 22  - Yi / exp(-1/2 * ((Xi - 7)/h)^2) - 1 + 16 - exp(-1/2 * ((Xi - 8)/h)^2) * 16  - Yi / exp(-1/2 * ((Xi - 8)/h)^2) - 1 + 10 - exp(-1/2 * ((Xi - 9)/h)^2) * 10  - Yi / exp(-1/2 * ((Xi - 9)/h)^2) - 1 + 18 - exp(-1/2 * ((Xi - 10)/h)^2) * 18  - Yi / exp(-1/2 * ((Xi - 10)/h)^2) - 1 + 26 - exp(-1/2 * ((Xi - 10)/h)^2) * 26  - Yi / exp(-1/2 * ((Xi - 10)/h)^2) - 1 + 34 - exp(-1/2 * ((Xi - 10)/h)^2) * 34  - Yi / exp(-1/2 * ((Xi - 10)/h)^2) - 1 + 17 - exp(-1/2 * ((Xi - 11)/h)^2) * 17  - Yi / exp(-1/2 * ((Xi - 11)/h)^2) - 1 + 28 - exp(-1/2 * ((Xi - 11)/h)^2) * 28  - Yi / exp(-1/2 * ((Xi - 11)/h)^2) - 1 + 14 - exp(-1/2 * ((Xi - 12)/h)^2) * 14  - Yi / exp(-1/2 * ((Xi - 12)/h)^2) - 1 + 20 - exp(-1/2 * ((Xi - 12)/h)^2) * 20  - Yi / exp(-1/2 * ((Xi - 12)/h)^2) - 1 + 24 - exp(-1/2 * ((Xi - 12)/h)^2) * 24  - Yi / exp(-1/2 * ((Xi - 12)/h)^2) - 1 + 28 - exp(-1/2 * ((Xi - 12)/h)^2) * 28  - Yi / exp(-1/2 * ((Xi - 12)/h)^2) - 1 + 26 - exp(-1/2 * ((Xi - 13)/h)^2) * 26  - Yi / exp(-1/2 * ((Xi - 13)/h)^2) - 1 + 34 - exp(-1/2 * ((Xi - 13)/h)^2) * 34  - Yi / exp(-1/2 * ((Xi - 13)/h)^2) - 1 + 34 - exp(-1/2 * ((Xi - 13)/h)^2) * 34  - Yi / exp(-1/2 * ((Xi - 13)/h)^2) - 1 + 46 - exp(-1/2 * ((Xi - 13)/h)^2) * 46  - Yi / exp(-1/2 * ((Xi - 13)/h)^2) - 1 + 26 - exp(-1/2 * ((Xi - 14)/h)^2) * 26  - Yi / exp(-1/2 * ((Xi - 14)/h)^2) - 1 + 36 - exp(-1/2 * ((Xi - 14)/h)^2) * 36  - Yi / exp(-1/2 * ((Xi - 14)/h)^2) - 1 + 60 - exp(-1/2 * ((Xi - 14)/h)^2) * 60  - Yi / exp(-1/2 * ((Xi - 14)/h)^2) - 1 + 80 - exp(-1/2 * ((Xi - 14)/h)^2) * 80  - Yi / exp(-1/2 * ((Xi - 14)/h)^2) - 1 + 20 - exp(-1/2 * ((Xi - 15)/h)^2) * 20  - Yi / exp(-1/2 * ((Xi - 15)/h)^2) - 1 + 26 - exp(-1/2 * ((Xi - 15)/h)^2) * 26  - Yi / exp(-1/2 * ((Xi - 15)/h)^2) - 1 + 54 - exp(-1/2 * ((Xi - 15)/h)^2) * 54  - Yi / exp(-1/2 * ((Xi - 15)/h)^2) - 1 + 32 - exp(-1/2 * ((Xi - 16)/h)^2) * 32  - Yi / exp(-1/2 * ((Xi - 16)/h)^2) - 1 + 40 - exp(-1/2 * ((Xi - 16)/h)^2) * 40  - Yi / exp(-1/2 * ((Xi - 16)/h)^2) - 1 + 32 - exp(-1/2 * ((Xi - 17)/h)^2) * 32  - Yi / exp(-1/2 * ((Xi - 17)/h)^2) - 1 + 40 - exp(-1/2 * ((Xi - 17)/h)^2) * 40  - Yi / exp(-1/2 * ((Xi - 17)/h)^2) - 1 + 50 - exp(-1/2 * ((Xi - 17)/h)^2) * 50  - Yi / exp(-1/2 * ((Xi - 17)/h)^2) - 1 + 42 - exp(-1/2 * ((Xi - 18)/h)^2) * 42  - Yi / exp(-1/2 * ((Xi - 18)/h)^2) - 1 + 56 - exp(-1/2 * ((Xi - 18)/h)^2) * 56  - Yi / exp(-1/2 * ((Xi - 18)/h)^2) - 1 + 76 - exp(-1/2 * ((Xi - 18)/h)^2) * 76  - Yi / exp(-1/2 * ((Xi - 18)/h)^2) - 1 + 84 - exp(-1/2 * ((Xi - 18)/h)^2) * 84  - Yi / exp(-1/2 * ((Xi - 18)/h)^2) - 1 + 36 - exp(-1/2 * ((Xi - 19)/h)^2) * 36  - Yi / exp(-1/2 * ((Xi - 19)/h)^2) - 1 + 46 - exp(-1/2 * ((Xi - 19)/h)^2) * 46  - Yi / exp(-1/2 * ((Xi - 19)/h)^2) - 1 + 68 - exp(-1/2 * ((Xi - 19)/h)^2) * 68  - Yi / exp(-1/2 * ((Xi - 19)/h)^2) - 1 + 32 - exp(-1/2 * ((Xi - 20)/h)^2) * 32  - Yi / exp(-1/2 * ((Xi - 20)/h)^2) - 1 + 48 - exp(-1/2 * ((Xi - 20)/h)^2) * 48  - Yi / exp(-1/2 * ((Xi - 20)/h)^2) - 1 + 52 - exp(-1/2 * ((Xi - 20)/h)^2) * 52  - Yi / exp(-1/2 * ((Xi - 20)/h)^2) - 1 + 56 - exp(-1/2 * ((Xi - 20)/h)^2) * 56  - Yi / exp(-1/2 * ((Xi - 20)/h)^2) - 1 + 64 - exp(-1/2 * ((Xi - 20)/h)^2) * 64  - Yi / exp(-1/2 * ((Xi - 20)/h)^2) - 1 + 66 - exp(-1/2 * ((Xi - 22)/h)^2) * 66  - Yi / exp(-1/2 * ((Xi - 22)/h)^2) - 1 + 54 - exp(-1/2 * ((Xi - 23)/h)^2) * 54  - Yi / exp(-1/2 * ((Xi - 23)/h)^2) - 1 + 70 - exp(-1/2 * ((Xi - 24)/h)^2) * 70  - Yi / exp(-1/2 * ((Xi - 24)/h)^2) - 1 + 92 - exp(-1/2 * ((Xi - 24)/h)^2) * 92  - Yi / exp(-1/2 * ((Xi - 24)/h)^2) - 1 + 93 - exp(-1/2 * ((Xi - 24)/h)^2) * 93  - Yi / exp(-1/2 * ((Xi - 24)/h)^2) - 1 + 120 - exp(-1/2 * ((Xi - 24)/h)^2) * 120  - Yi / exp(-1/2 * ((Xi - 24)/h)^2) - 1 + 85 - exp(-1/2 * ((Xi - 25)/h)^2) * 85  - Yi / exp(-1/2 * ((Xi - 25)/h)^2) - 1 )^2"

なんか幾何学模様に見えてきたぞ。。。