Rで2式の交点を求める

Question

#Rで求めた２つの回帰曲線との交点座標を求めたいです。

以下のようなcsvファイルデータから左と回答する確率(ans_left)と、右と回答する確率(ans_right)を、degを説明変数とするロジスティック回帰を求めました。
そこでそれら回帰曲線らの交点を求めたいです。

	imagename	ans_left	ans_direct	ans_right	direction	deg
0	any_0_+20	0	0	1	0	20
1	any_0_-15	1	0	0	0	-15
2	any_+30_-5	0	1	0	30	-5
...

• ans_left：左と回答したかどうか(0 or 1)←従属変数
• ans_right：右と回答したかどうか(0 or 1)←従属変数
• deg： 角度(-20~20) ←説明変数

• 赤曲線：glm()で求めたロジスティック回帰曲線1
• 緑曲線：glm()で求めたロジスティック回帰曲線2
• 青曲線: 1- 赤曲線 - 緑曲線

求めたいこと

• 赤曲線と青曲線の交点座標
• 緑曲線と青曲線の交点座標
• 青曲線の最大値とその座標

現在のコード

曲線をプロットすることは出来ましたが、交点座標の求め方がわかりません。

HEAD_L = -30
HEAD_D = 0
HEAD_R = 30

# plotする図のx軸(deg)の上下限 #
xMin = -30 # x軸の下限
xMax = +30 # x軸の上限


# データの読み込み #
df_long <- read.csv('log/gomi.csv',header=T) # 1列目(header)=列名(True)
df_short <- read.csv('log/0826_2/0826_2_short.csv',header=T)

# 関数宣言 #
f_ansL <- function(ansData){ #ansData:列名"direction"の値=-30の行の、列3("ans_left"), 7("deg")のデータ
  ans <- glm(ans_left~deg, data = ansData, family=binomial(logit))
  summary(ans) # 要約
  return (ans)
}

f_ansR <- function(ansData){ #ansData:列名"direction"の値=-30の行の、列3("ans_left"), 7("deg")のデータ
  ans <- glm(ans_right~deg, data = ansData, family=binomial(logit))
  summary(ans) # 要約
  return (ans)
}

long_headL_ansL_data <- df_long[df_long$direction== HEAD_L, c(3,7)] # 列名"direction"の値=-30の行の、列3("ans_left"), 7("deg")の値を取得 # 列の値が条件を満たす行データだけ抽出
long_headL_ansL = f_ansL(long_headL_ansL_data)

long_headL_ansR_data <- df_long[df_long$direction== HEAD_L, c(5,7)] # 列名"direction"の値=-30の行の、列3("ans_left"), 7("deg")の値を取得
long_headL_ansR = f_ansR(long_headL_ansR_data)



y <- tapply(long_headL_ansL_data$ans_left, long_headL_ansL_data$deg, mean) 
x <- names(y)
points(x, y,col=2)

par(new = TRUE) # 図の重ね合わせを許可

y_long_headL_ansR <- predict(long_headL_ansR, data.frame(deg), type = "response") 
plot(deg, y_long_headL_ansR, type = "l", xlim = c(xMin, xMax), ylim = c(0, 1), ann = F , col = 3)


deg <- seq(xMin, xMax, length = 100)

y_long_headL_ansL <- predict(long_headL_ansL, data.frame(deg), type = "response") 
png("logisticGlm0000.png")
plot(deg, y_long_headL_ansL, type = "l", xlim = c(xMin, xMax), ylim = c(0, 1), xlab="direcion of gaze[deg]", ylab="Prop. of resposes", col = 2)


y <- tapply(long_headL_ansR_data$ans_right, long_headL_ansR_data$deg, mean)
x <- names(y)
points(x, y,col=3)
par(new = TRUE) # 図の重ね合わせを許可
plot(deg, 1- y_long_headL_ansL - y_long_headL_ansR, type = "l", xlim = c(xMin, xMax), ylim = c(0, 1), ann = F , col = 4)

Answer 1

結局Rでの解き方がイマイチわからず...
Pythonに連立方程式の解を求めるfsolve関数というものがあり、リファレンスや記事も充実していたことから、Rで回帰式導出→pythonで交点導出で一旦解決しました。
(これならいっそ回帰式導出もPythonで再実装すべきかなと思ったり.....)

11.scipyの基本と応用
こちらのサイトのexampele-02を参考にさせて頂きました。
↑こちらでは交点が2つある場合の求め方と、グラフへのプロットまで丁寧に記述されていました！

私が求たいのは一点の交点だけだったので、必要な部分だけ真似してコーディングさせて頂きました。

式1と式2の1交点を求める

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve



x = np.linspace(-20, 20, 1000)  # -20~20に10000の要素を持つ等差数列を生成
def func(p):
    x, y = p[0], p[1]
    f = np.zeros(2)  # 0を初期値とする2要素の配列を生成
    f[0] = 式1  # y = 3x + 1　のときは 3*x + 1 - y と記述
    f[1] = 式2  
    return f

o1 = fsolve(func, [1, 1])
o1[0] # 交点のx座標
o1[1] # 交点のy座標

ちなみに
Pythonで方程式を解く方法（SciPy、ニュートン法、二分法による計算）
によると方程式を解く方法も複数あるようです。
ニュートン法や二分法の方が解導出が早いようですが、今回は大した計算ではないし、処理速度は気にしていなかったので実装が一番手軽なSciPyによる計算を行いました。

Answer 2

解析的に求めたいのでしょうか？数値的に求めたいのでしょうか
1/(1+exp(-ax_1+b)) = 1/(1+exp(-cx_2+d)なので解析的にも求められるかと思います

数値的にでればx-yの座標がありますので、ループしてyの値がイコールか微小な一定範囲内となるxの値を出せばよいかと思います。

Comments

[funyo]

回答ありがとうございます。
導出された回帰式をもとにrプログラムで解析的に求める場合は、関数nleqslvを用いて求めるということでしょうか？
(rの関数に関して知識が浅いため稚拙な質問で申し訳ありません。間違いや他に推奨される関数があれば教えていただきたいです。)

[aokikenichi]

いえ、プログラムでなくて手計算でもよろしいかと。
大量にあるとのことでしょうか？

[funyo]

確かに手計算でも可能ですね。
ただ、いくつか同じ解析を行うのでプログラム内で自動化出来たらベストです。

[aokikenichi]

私は、nleqslvほかRで数式解析をしたことがないので、申し訳ないですが分かりません。

私であれば数値で行きます。