ひゃー。
> 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"
なんか幾何学模様に見えてきたぞ。。。