Md 根を求める・近似曲線を引く

提供: StatGenetKyotoU
移動: 案内検索
根を求める・近似曲線を引く
カテゴリ 名称 説明
医学生物学トピック データを読む,用量反応曲線
数学トピック 最適化、最尤法、最小二乗法 もっともあてはまりの良い値の近似値を出す
R関数・Rスキル ソースを読む、uniroot(),optim(),glm() 収束、近似

t/sin(t)=a ({\frac  {t}{\sin(t)}}=a)を解く

    • 双曲幾何的問題として、ひだの図示のために上記の等式を解きたい(こちら)
L<-1
a<-1

fx<-function(x,a){
	x/sin(x)-a
}
interval<-c(0,pi/2)
uniroot(fx,interval=c(10^(-10),pi/2),a=1.1)

X<-1
Y<-1

Ny<-100
Nx<-100

ys<-seq(from=0,to=1,length=Ny)*Y
Xs<-Ys<-Zs<-c()
Ls<-a*ys+L

for(i in 1:Ny){
	L2<-Ls[i]
	if(L2==L){
		Xs<-c(Xs,seq(from=0,to=1,length=Nx)*X)
		Ys<-c(Ys,rep(ys[i],Nx))
		Zs<-c(Zs,rep(0,Nx))
	}else{
		t<-uniroot(fx,interval=c(10^(-10),pi/2),a=L2/L)[[1]]
		ts<-seq(from=-t/2,to=t/2,length=Nx)
		r<-L/(2*sin(t/2))
		Xc<-L/2
		Zc<-L/2/tan(t/2)
		Ys<-c(Ys,rep(ys[i],Nx))
		Xs<-c(Xs,r*sin(ts)+Xc)
		Zs<-c(Zs,r*cos(ts)-Zc)
	}
}
xlim<-ylim<-zlim<-range(c(Xs,Ys,Zs))
plot3d(Xs,Ys,Zs,xlim=xlim,ylim=ylim,zlim=zlim)

データによくあてはまる多項式を推定する

# Function to estimate polynomial coefficients and coordinates estimated
# 多項式近似の係数と、推定座標を返す関数
glmPolynomial<-function(y,x,dg){# y:dependent values (従属値列),x:descriptive values (説明変数),dg:degrees of polynomials (多項式の次数)
 # Xs: Matrix consisted of x^1,x^2,...x^{dg}
 # Xs: xの1,2,...,dg 乗の値でできた行列
 Xs<-matrix(0,length(x),dg) 
 for(i in 1:dg){
  Xs[,i]<-x^i
 }
 # glm() is generalized linear regression function in R
 # glm() はRの一般化線形回帰関数
 fit<-glm(y~Xs[,1:dg]) 
 # beta: coefficients estimated
 # beta: 推定係数ベクトル
 beta<-fit$coefficients
 # Xs2 has an additional column for x^0
 # Xs2 には、x^0のための列を追加
 Xs2<-cbind(rep(1,length(x)),Xs)
 # predY is a matrix of values for individual degrees of x
 # predY は xの各次数の項
 predY<-t(t(Xs2)*beta)
 # apply(): Sum up all columns corresponding to 0 to dg degrees
 # apply()関数で0:dg次数のすべてを足し合わせる
 sumupPredY<-apply(predY,1,sum)
 list(beta=beta,x=x,predY=sumupPredY)
}

# Generate data
x<-sort(runif(100))
y<-runif(100)+x^3*2+x^4*5

xlim<-c(min(x),max(x))
ylim<-c(min(y),max(y))
plot(x,y,xlim=xlim,ylim=ylim,col="red")

# GLM is to be applied to multiple degrees (dgs)
# GLMを複数の次数に適用
dgs<-1:5
for(i in dgs){
 outglm<-glmPolynomial(y,x,i)
 par(new=TRUE)
 # gray() changes contrast 
 # gray() 関数は線の濃淡を変える
 plot(outglm$x,outglm$predY,xlim=xlim,ylim=ylim,col=gray(i/(max(dgs)*1.5)),type="l")
}

用量反応曲線を引く

Rの関数の中身を読む

その他の記事