hack のためのネタ帳, etc,,,

サンプルデータ

ここでは y = x^2 についてノイズ入りのデータを生成してみる。
set.seed(0)
d = list()
d$x =  1:10   ; d$x = d$x + rnorm(length(d$x))
d$y = (1:10)^2; d$y = d$y + rnorm(length(d$y))
> d
$x
 [1]  2.262954  1.673767  4.329799  5.272429  5.414641  4.460050  6.071433  7.705280  8.994233 12.404653

$y
 [1]  1.763593  3.200991  7.852343 15.710538 24.700785 35.588489 49.252223 63.108079 81.435683 98.762462

曲線の描画

近似式の左辺値を推定するために必要な右辺値を適当な密度で生成して predict() に与えることで近似式に当てはめた左辺値を生成し、それらを lines() 等に与えると良いらしい。
詳細は、以後のサンプルを参照。

lm() を使う方法

lm() は Fitting Linear Models つまり線形モデルによる近似を行うらしい。
formula に近似式、data にデータを与える。
近似式の与え方は ~ の仕組みが良く分かってないのだが、与えられた data の項目を ~ で左辺と右辺に分けて与えると、右辺の N 次式に対して左辺の値を求めるために適当な係数が得られる。
lm() の場合、近似式には係数を与える必要はないようだが、高次の項は I(x^2) のようにラップしておかないとなぜかすべて一次の項として扱われ一次の項の係数しか得られない。
また、ゼロ次の項が必ず入ってしまうようで、これをネグる方法が分からない。
この辺りの話は、どこかに解説があるはずなので要捜索。
あとは lm() の結果と右辺で用いたパラメータに割り当てる値を predict() に与えると、それらを近似式に当てはめた際の左辺値が得られる。
以下の例では、d$x の最小値から最大値までの間に 100 点生成して lm() の結果と共に predict() へ与えることで、近似曲線を描くのに十分な数の点列を得て lines() で描画している。
data に与えるのは data.frame() の方が良い気もするのだが、回帰分析の結果も一緒に放り込んでおきたいので、list() にしてパラメータ名を揃えて放り込むことで手間を省いている。
d.lm = list()
d.lm$obj = lm(formula = y ~ I(x^2), data = d)
d.lm$x = seq(min(d$x), max(d$x), length=100)
d.lm$y = predict(d.lm$obj, d.lm)
> d.lm
$obj

Call:
lm(formula = y ~ I(x^2), data = d)

Coefficients:
(Intercept)       I(x^2)  
     8.2542       0.6875  


$x
  [1]  1.673767  1.782159  1.890552  1.998945  2.107338  2.215731  2.324123  2.432516  2.540909  2.649302  2.757695  2.866087  2.974480  3.082873
 [15]  3.191266  3.299659  3.408051  3.516444  3.624837  3.733230  3.841623  3.950015  4.058408  4.166801  4.275194  4.383587  4.491979  4.600372
 [29]  4.708765  4.817158  4.925551  5.033943  5.142336  5.250729  5.359122  5.467514  5.575907  5.684300  5.792693  5.901086  6.009478  6.117871
 [43]  6.226264  6.334657  6.443050  6.551442  6.659835  6.768228  6.876621  6.985014  7.093406  7.201799  7.310192  7.418585  7.526978  7.635370
 [57]  7.743763  7.852156  7.960549  8.068942  8.177334  8.285727  8.394120  8.502513  8.610906  8.719298  8.827691  8.936084  9.044477  9.152870
 [71]  9.261262  9.369655  9.478048  9.586441  9.694834  9.803226  9.911619 10.020012 10.128405 10.236797 10.345190 10.453583 10.561976 10.670369
 [85] 10.778761 10.887154 10.995547 11.103940 11.212333 11.320725 11.429118 11.537511 11.645904 11.754297 11.862689 11.971082 12.079475 12.187868
 [99] 12.296261 12.404653

$y
        1         2         3         4         5         6         7         8         9        10        11        12        13        14        15 
 10.18016  10.43768  10.71136  11.00120  11.30718  11.62933  11.96762  12.32207  12.69268  13.07944  13.48235  13.90142  14.33664  14.78802  15.25555 
       16        17        18        19        20        21        22        23        24        25        26        27        28        29        30 
 15.73923  16.23907  16.75506  17.28721  17.83551  18.39997  18.98058  19.57734  20.19026  20.81933  21.46456  22.12594  22.80348  23.49717  24.20701 
       31        32        33        34        35        36        37        38        39        40        41        42        43        44        45 
 24.93301  25.67516  26.43347  27.20793  27.99854  28.80531  29.62823  30.46731  31.32254  32.19393  33.08147  33.98516  34.90501  35.84102  36.79317 
       46        47        48        49        50        51        52        53        54        55        56        57        58        59        60 
 37.76148  38.74595  39.74657  40.76334  41.79627  42.84535  43.91059  44.99198  46.08953  47.20323  48.33308  49.47909  50.64125  51.81957  53.01404 
       61        62        63        64        65        66        67        68        69        70        71        72        73        74        75 
 54.22466  55.45144  56.69437  57.95346  59.22870  60.52010  61.82765  63.15135  64.49121  65.84722  67.21939  68.60771  70.01219  71.43282  72.86960 
       76        77        78        79        80        81        82        83        84        85        86        87        88        89        90 
 74.32254  75.79163  77.27688  78.77828  80.29583  81.82954  83.37940  84.94542  86.52759  88.12592  89.74040  91.37103  93.01782  94.68077  96.35986 
       91        92        93        94        95        96        97        98        99       100 
 98.05512  99.76652 101.49408 103.23779 104.99766 106.77369 108.56586 110.37419 112.19868 114.03932 

plot(d)
lines(d.lm)

glm() を使う方法

glm() は Fitting Generalized Linear Modesls つまり、一般化線形モデルによる近似を行うらしい。
基本的には lm() と同じだが、start=c(1,1) みたいにして、係数の初期値を与えることが出来たり、他にも一般化とついてるだけあっていろいろと設定ができるようだ。
このくらい単純な例だとあまり有難味はなさそう。
d.glm = list()
d.glm$obj = glm(formula = y ~ I(x^2), data = d)
d.glm$x = seq(min(d$x), max(d$x), length=100)
d.glm$y = predict(d.glm$obj, d.glm)
> d.glm
$obj

Call:  glm(formula = y ~ I(x^2), data = d)

Coefficients:
(Intercept)       I(x^2)  
     8.2542       0.6875  

Degrees of Freedom: 9 Total (i.e. Null);  8 Residual
Null Deviance:	    10450 
Residual Deviance: 1646 	AIC: 85.41

$x
  [1]  1.673767  1.782159  1.890552  1.998945  2.107338  2.215731  2.324123  2.432516  2.540909  2.649302  2.757695  2.866087  2.974480  3.082873
 [15]  3.191266  3.299659  3.408051  3.516444  3.624837  3.733230  3.841623  3.950015  4.058408  4.166801  4.275194  4.383587  4.491979  4.600372
 [29]  4.708765  4.817158  4.925551  5.033943  5.142336  5.250729  5.359122  5.467514  5.575907  5.684300  5.792693  5.901086  6.009478  6.117871
 [43]  6.226264  6.334657  6.443050  6.551442  6.659835  6.768228  6.876621  6.985014  7.093406  7.201799  7.310192  7.418585  7.526978  7.635370
 [57]  7.743763  7.852156  7.960549  8.068942  8.177334  8.285727  8.394120  8.502513  8.610906  8.719298  8.827691  8.936084  9.044477  9.152870
 [71]  9.261262  9.369655  9.478048  9.586441  9.694834  9.803226  9.911619 10.020012 10.128405 10.236797 10.345190 10.453583 10.561976 10.670369
 [85] 10.778761 10.887154 10.995547 11.103940 11.212333 11.320725 11.429118 11.537511 11.645904 11.754297 11.862689 11.971082 12.079475 12.187868
 [99] 12.296261 12.404653

$y
        1         2         3         4         5         6         7         8         9        10        11        12        13        14        15 
 10.18016  10.43768  10.71136  11.00120  11.30718  11.62933  11.96762  12.32207  12.69268  13.07944  13.48235  13.90142  14.33664  14.78802  15.25555 
       16        17        18        19        20        21        22        23        24        25        26        27        28        29        30 
 15.73923  16.23907  16.75506  17.28721  17.83551  18.39997  18.98058  19.57734  20.19026  20.81933  21.46456  22.12594  22.80348  23.49717  24.20701 
       31        32        33        34        35        36        37        38        39        40        41        42        43        44        45 
 24.93301  25.67516  26.43347  27.20793  27.99854  28.80531  29.62823  30.46731  31.32254  32.19393  33.08147  33.98516  34.90501  35.84102  36.79317 
       46        47        48        49        50        51        52        53        54        55        56        57        58        59        60 
 37.76148  38.74595  39.74657  40.76334  41.79627  42.84535  43.91059  44.99198  46.08953  47.20323  48.33308  49.47909  50.64125  51.81957  53.01404 
       61        62        63        64        65        66        67        68        69        70        71        72        73        74        75 
 54.22466  55.45144  56.69437  57.95346  59.22870  60.52010  61.82765  63.15135  64.49121  65.84722  67.21939  68.60771  70.01219  71.43282  72.86960 
       76        77        78        79        80        81        82        83        84        85        86        87        88        89        90 
 74.32254  75.79163  77.27688  78.77828  80.29583  81.82954  83.37940  84.94542  86.52759  88.12592  89.74040  91.37103  93.01782  94.68077  96.35986 
       91        92        93        94        95        96        97        98        99       100 
 98.05512  99.76652 101.49408 103.23779 104.99766 106.77369 108.56586 110.37419 112.19868 114.03932 

plot(d)
lines(d.glm)

nls() を使う方法

nls() は Nonelinear Least Square つまり非線形最小二乗法を行うらしい。
これは、lm() や glm() と違って formula には近似式を係数込みで与え、start で係数の初期値を与える必要があるが、そのおかげでゼロ次の項をネグることが出来る。
d.nls = list()
d.nls$obj = nls(formula = y ~ a*x^2, data = d, start = c(a=1))
d.nls$x = seq(min(d$x), max(d$x), length=100)
d.nls$y = predict(d.nls$obj, d.nls)
> d.nls
$obj
Nonlinear regression model
  model: y ~ a * x^2
   data: d
     a 
0.7831 
 residual sum-of-squares: 1984

Number of iterations to convergence: 1 
Achieved convergence tolerance: 3.25e-09

$x
  [1]  1.673767  1.782159  1.890552  1.998945  2.107338  2.215731  2.324123  2.432516  2.540909  2.649302  2.757695  2.866087  2.974480  3.082873
 [15]  3.191266  3.299659  3.408051  3.516444  3.624837  3.733230  3.841623  3.950015  4.058408  4.166801  4.275194  4.383587  4.491979  4.600372
 [29]  4.708765  4.817158  4.925551  5.033943  5.142336  5.250729  5.359122  5.467514  5.575907  5.684300  5.792693  5.901086  6.009478  6.117871
 [43]  6.226264  6.334657  6.443050  6.551442  6.659835  6.768228  6.876621  6.985014  7.093406  7.201799  7.310192  7.418585  7.526978  7.635370
 [57]  7.743763  7.852156  7.960549  8.068942  8.177334  8.285727  8.394120  8.502513  8.610906  8.719298  8.827691  8.936084  9.044477  9.152870
 [71]  9.261262  9.369655  9.478048  9.586441  9.694834  9.803226  9.911619 10.020012 10.128405 10.236797 10.345190 10.453583 10.561976 10.670369
 [85] 10.778761 10.887154 10.995547 11.103940 11.212333 11.320725 11.429118 11.537511 11.645904 11.754297 11.862689 11.971082 12.079475 12.187868
 [99] 12.296261 12.404653

$y
  [1]   2.193858   2.487206   2.798956   3.129107   3.477659   3.844613   4.229968   4.633725   5.055882   5.496441   5.955402   6.432764   6.928527
 [14]   7.442691   7.975257   8.526224   9.095593   9.683362  10.289533  10.914106  11.557080  12.218455  12.898232  13.596410  14.312989  15.047969
 [27]  15.801351  16.573135  17.363319  18.171905  18.998892  19.844281  20.708071  21.590262  22.490855  23.409849  24.347245  25.303041  26.277239
 [40]  27.269839  28.280839  29.310242  30.358045  31.424250  32.508856  33.611863  34.733272  35.873082  37.031294  38.207907  39.402921  40.616337
 [53]  41.848153  43.098372  44.366991  45.654012  46.959434  48.283258  49.625483  50.986109  52.365137  53.762566  55.178396  56.612628  58.065261
 [66]  59.536295  61.025731  62.533568  64.059807  65.604446  67.167487  68.748930  70.348774  71.967019  73.603665  75.258713  76.932162  78.624013
 [79]  80.334265  82.062918  83.809972  85.575428  87.359285  89.161544  90.982204  92.821265  94.678728  96.554592  98.448857 100.361524 102.292592
 [92] 104.242061 106.209932 108.196204 110.200877 112.223952 114.265428 116.325305 118.403584 120.500264

plot(d)
lines(d.nls)

参考

コメントをかく


「http://」を含む投稿は禁止されています。

利用規約をご確認のうえご記入下さい

Wiki内検索

フリーエリア

編集にはIDが必要です