{{category R}} R !!!非線形回帰分析 {{outline}} ---- !!nls() nls(formula,start,trace) nls(式, パラメタの初期値) !残差のプロット plot(resid(モデル)) !References http://cse.naro.affrc.go.jp/takezawa/r-tips/r/73.html !!例:Syntactic Tree Fragmentsがzipfの法則にしたがうか {{pre fragments.dat <- read.table("fragments.dat.txt") head(fragments.dat) ## fragment rank freq ## 1 (PP (IN ) (NP )) 1 6257 ## 2 (NP (PRP )) 2 4242 ## 3 (S (NP ) (VP )) 3 4079 ## 4 (S (VP )) 4 3714 ## 5 (DT the) 5 3450 ## 6 (NP (DT ) (NN )) 6 3260 frag.short <- head(fragments.dat, 500) dim(frag.short) ## [1] 500 3 str(frag.short) ## 'data.frame': 500 obs. of 3 variables: ## $ fragment: chr "(PP (IN ) (NP ))" "(NP (PRP ))" "(S (NP ) (VP ))" "(S (VP ))" ... ## $ rank : int 1 2 3 4 5 6 7 8 9 10 ... ## $ freq : int 6257 4242 4079 3714 3450 3260 3101 3085 2850 2569 ... plot(frag.short$rank, frag.short$freq, xlim=c(0,500), ylim=c(0,7000)) }} {{ref_image rank.freq.png}} !Zipfの法則のモデル式 {{ref_image zipf.png}} {{pre fit <- nls(freq ~ a/(rank^b), start=list(a=6257,b=1), data=frag.short, trace=T) 1.110042e+08 (4.97e+00): par = (6257 1) 4.156852e+07 (1.59e+00): par = (6095.586 0.6873339) 1.741337e+07 (5.77e-01): par = (6665.465 0.6063987) 1.308021e+07 (2.59e-02): par = (7571.729 0.5972984) 1.307277e+07 (4.32e-03): par = (7599.867 0.6002608) 1.307257e+07 (8.20e-04): par = (7590.96 0.5997113) 1.307256e+07 (1.55e-04): par = (7592.631 0.5998156) 1.307256e+07 (2.94e-05): par = (7592.314 0.5997958) 1.307256e+07 (5.58e-06): par = (7592.374 0.5997996) Formula: freq ~ a/(rank^b) Parameters: Estimate Std. Error t value Pr(>|t|) a 7.592e+03 1.210e+02 62.74 <2e-16 *** b 5.998e-01 5.694e-03 105.34 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 162 on 498 degrees of freedom Number of iterations to convergence: 8 Achieved convergence tolerance: 5.576e-06 plot(frag.short$rank, frag.short$freq, xlim=c(0,500), ylim=c(0,7000)) lines(frag.short$rank, fitted(fit), col="red", lwd=3) }} {{ref_image zipf.fit.png}} !得られたパラメタを使ったモデル式 {{ref_image zipf.tf.png}} !どのくらいフィットしているか、残差分析プロット plot(resid(fit)) {{ref_image zipf.resid.png}}