ベストフィットの多項式とプロットのドロップダウンラインを

Question

私はWindows 10でR 3.3.1（64ビット）を使用しています。私は2次多項式に適合するx-yデータセットを持っています。私は、y = 4でxに最もよくフィットする多項式を解き、y = 4からx軸にドロップダウンラインをプロットしたいと思います。

これは、データフレームV1内のデータを生成します。

v1 <- structure(list(x = c(-5.2549, -3.4893, -3.5909, -2.5546, -3.7247, 
-5.1733, -3.3451, -2.8993, -2.6835, -3.9495, -4.9649, -2.8438, 
-4.6926, -3.4768, -3.1221, -4.8175, -4.5641, -3.549, -3.08, -2.4153, 
-2.9882, -3.4045, -4.6394, -3.3404, -2.6728, -3.3517, -2.6098, 
-3.7733, -4.051, -2.9385, -4.5024, -4.59, -4.5617, -4.0658, -2.4986, 
-3.7559, -4.245, -4.8045, -4.6615, -4.0696, -4.6638, -4.6505, 
-3.7978, -4.5649, -5.7669, -4.519, -3.8561, -3.779, -3.0549, 
-3.1241, -2.1423, -3.2759, -4.224, -4.028, -3.3412, -2.8832, 
-3.3866, -0.1852, -3.3763, -4.317, -5.3607, -3.3398, -1.9087, 
-4.431, -3.7535, -3.2545, -0.806, -3.1419, -3.7269, -3.4853, 
-4.3129, -2.8891, -3.0572, -5.3309, -2.5837, -4.1128, -4.6631, 
-3.4695, -4.1045, -7.064, -5.1681, -6.4866, -2.7522, -4.6305, 
-4.2957, -3.7552, -4.9482, -5.6452, -6.0302, -5.3244, -3.9819, 
-3.8123, -5.3085, -5.6096, -6.4557), y = c(0.99, 0.56, 0.43, 
2.31, 0.31, 0.59, 0.62, 1.65, 2.12, 0.1, 0.24, 1.68, 0.09, 0.59, 
1.23, 0.4, 0.36, 0.49, 1.41, 3.29, 1.22, 0.56, 0.1, 0.67, 2.38, 
0.43, 1.56, 0.07, 0.08, 1.53, -0.01, 0.12, 0.1, 0.04, 3.42, 0.23, 
0, 0.34, 0.15, 0.03, 0.19, 0.17, 0.2, 0.09, 2.3, 0.07, 0.15, 
0.18, 1.07, 1.21, 3.4, 0.8, -0.04, 0.02, 0.74, 1.59, 0.71, 10.64, 
0.64, -0.01, 1.06, 0.81, 4.58, 0.01, 0.14, 0.59, 7.35, 0.63, 
0.17, 0.38, -0.08, 1.1, 0.89, 0.94, 1.52, 0.01, 0.1, 0.38, 0.02, 
7.76, 0.72, 4.1, 1.36, 0.13, -0.02, 0.13, 0.42, 1.49, 2.64, 1.01, 
0.08, 0.22, 1.01, 1.53, 4.39)), .Names = c("x", "y"), class = "data.frame", row.names = c(NA, 
-95L)) 

はここで、x対yのプロットベストフィット多項式をプロットし、Y = 4で線を描画するためのコードです。

> attach(v1) 
> # simple x-y plot of the data 
> plot(x,y, pch=16) 
> # 2nd order polynomial fit 
> fit2 <- lm(y~poly(x,2,raw=TRUE)) 
> summary(fit2) 
> # generate range of numbers for plotting polynomial 
> xx <- seq(-8,0, length=50) 
> # overlay best fit polynomial 
>lines(xx, predict(fit2, data.frame(x=xx)), col="blue") 
> # add horizontal line at y=4 
> abline(h=4, col="red") 
> 

これは、Yの周りのX -2、-6.5で4 =プロットから明らかだが、私は実際にこれらの値についての回帰多項式を解決したいです。

理想的には、赤青の線の交点からx軸（つまり、2つのy = 4の解で終了するプロットの垂直線）にドロップする線が理想的です。それが不可能なのであれば、適切なx解の値をとっている限り、プロットの上にある古い古い垂直線に満足しています。

このグラフは、y> 4のときに規格外となる部分を表しています。したがって、スペック内の部分を生成するx値の範囲をドロップダウンラインで強調したいと考えています。

Answer 1

あなたが値を計算するために次の式を使用することができます。

betas <- coef(fit2) # get coefficients 
betas[1] <- betas[1] - 4 # adjust intercept to look for values where y = 4 

# note degree increases, so betas[1] is c, etc. 
betas 
##    (Intercept) poly(x, 2, raw = TRUE)1 poly(x, 2, raw = TRUE)2 
##    8.7555833    6.0807302    0.7319848 

solns <- c((-betas[2] + sqrt(betas[2]^2 - 4 * betas[3] * betas[1]))/(2 * betas[3]), 
      (-betas[2] - sqrt(betas[2]^2 - 4 * betas[3] * betas[1]))/(2 * betas[3])) 

solns 
## poly(x, 2, raw = TRUE)1 poly(x, 2, raw = TRUE)1 
##    -1.853398    -6.453783 

segments(solns, -1, solns, 4, col = 'green') # add segments to graph 

（あなたがそれを見つけることができるかどうか）

はるかに簡単ですpolyroot：

polyroot(betas) 
## [1] -1.853398+0i -6.453783+0i 

それは複雑なベクトルを返すので、あなたはsegmentsにそれを渡したい場合はas.numericでそれをラップする必要があります。

Answer 2

あなたはあなただけの解析的にこの問題を解決するために、二次式を使用することができ、二次方程式

0.73198 * x^2 + 6.08073 * x + 12.75558 = 4 
OR 
0.73198 * x^2 + 6.08073 * x + 8.75558 = 0 

を持っています。 Rは、二つの根を与える：

(-6.08073 + sqrt(6.08073^2 -4*0.73198 * 8.75558))/(2 * 0.73198) 
[1] -1.853392 
(-6.08073 - sqrt(6.08073^2 -4*0.73198 * 8.75558))/(2 * 0.73198) 
[1] -6.453843 

abline（のV = Cの（-1.853392、-6.453843））

Answer 3

この単純2次多項式の解析解があることは絶対に理解できます。私が数値解を示す理由は、回帰の設定でこの質問をすることです。より複雑な回帰曲線を持つとき、数値解法は常に一般的にあなたの解法になります。以下では

私はuniroot機能を使用します。あなたがそれに慣れていないなら、最初に短い答え、Uniroot solution in Rを読んでください。

---

これはあなたのコードを生成プロットです。あなたはほとんどそこにいます。これは根の発見の問題です。unirootを数値で使用することができます。図から

f <- function (x) { 
    ## subtract 4 
    predict(fit2, newdata = data.frame(x = x)) - 4 
    } 

、2つの根、[-7, -6]内部1、[-3, -1]内の他があることは明らかである：の関数を定義してみましょう。我々は両方を見つけるためにunirootを使用します。

x1 <- uniroot(f, c(-7, -6))$root 
#[1] -6.453769 

x2 <- uniroot(f, c(-3, -1))$root 
#[1] -1.853406 

今、あなたは、x軸に至るまで、これらの点から垂直線をドロップすることができます。ここでは

y1 <- f(x1) + 4 ## add 4 back 
y2 <- f(x2) + 4 

abline(h = 0, col = 4) ## x-axis 
segments(x1, 0, x1, y1, lty = 2) 
segments(x2, 0, x2, y2, lty = 2) 

Answer 4

は、1つの以上のソリューションです、に基づいthis

attach(v1) 
fit2 = lm(y~poly(x,2,raw=TRUE)) 
xx = seq(-8,0, length=50) 

vector1 = predict(fit2, data.frame(x=xx)) 
vector2= replicate(length(vector1),4) 

# Find points where vector1 is above vector2. 
above = vector1 > vector2 

# Points always intersect when above=TRUE, then FALSE or reverse 
intersect.points = which(diff(above)!=0)  

# Find the slopes for each line segment. 
vector1.slopes = vector1[intersect.points+1] - vector1[intersect.points] 
vector2.slopes = vector2[intersect.points+1] - vector2[intersect.points] 

# Find the intersection for each segment. 
x.points = intersect.points + ((vector2[intersect.points] - vector1[intersect.points])/(vector1.slopes-vector2.slopes)) 
y.points = vector1[intersect.points] + (vector1.slopes*(x.points-intersect.points)) 

#Scale x.points to the axis value of xx 
x.points = xx[1] + ((x.points - 1)/(49))*(xx[50]-xx[1]) 

plot(xx, y = vector1, type= "l", col = "blue") 
points(x,y,pch = 20) 
lines(x = c(x.points[1],x.points[1]), y = c(0,y.points[1]), col='red') 
lines(x = c(x.points[2],x.points[2]), y = c(0,y.points[2]), col='red') 

Answer 5

すでに多くの解決策が提案されていますが、ここではもう1つです。

明らかに、多項式（二次）式a_0 + a_1.x + a_2.x^2 = 4を満たすxの値を見つけることに興味があります。ここで、a_0, a_1, a_2は適合多項式の係数です。我々は、標準的な二次方程式ax^2+bx+c=0として方程式を書き換えると次のように多項式回帰を取り付けた多項式の係数を用いてSridhar's式を用いて根を見つけることができる：

我々は、いくつかの数値的方法を使用することができ

a <- fit2$coefficients[3] 
b <- fit2$coefficients[2] 
c <- fit2$coefficients[1] - 4 

as.numeric((-b + sqrt(b^2-4*a*c))/(2*a)) 
#[1] -1.853398 
as.numeric((-b-+ sqrt(b^2-4*a*c))/(2*a)) 
#[1] -6.453783 

数値計算と理論解が一致しているように、Newton-Raphsonとして根を見つけることができます（高速数値計算方法がありますが、これは私たちの目的を解決し、かなり速いです。私のマシンで~160 msを使用します）。 。

a <- fit2$coefficients # fitted quadratic polynomial coefficients 

f <- function(x) { 
    as.numeric(a[1] + a[2]*x + a[3]*x^2-4) 
} 

df <- function(x) { 
    as.numeric(a[2] + 2*a[3]*x) 
} 

Newton.Raphson <- function(x0) { 
    eps <- 1e-6 
    x <- x0 
    while(TRUE) { 
    x <- x0 - f(x0)/df(x0) 
    if (abs(x - x0) < eps) { 
     return(x0) 
    } 
    x0 <- x 
    } 
} 

t1 <- Sys.time() 
x1 <- Newton.Raphson(-10) 
x2 <- Newton.Raphson(10) 
x1 
#[1] -6.453783 
x2 
#[1] -1.853398 
s2 
print(paste('time taken to compute the roots:' ,Sys.time() - t1)) 
#[1] "time taken to compute the roots: 0.0160109996795654" 
points(x1, 4, pch=19, col='green') 
points(x2, 4, pch=19, col='green') 
abline(v=x1, col='green') 
abline(v=x2, col='green')