冪乗近似について 冪乗近似の傾向にあるデータを近似直線でY=aX^bの形にしようと思ったのですが l

冪乗近似について　
冪乗近似の傾向にあるデータを近似直線でY=aX^bの形にしようと思ったのですが　lnx-lny のグラフをとったのち計算すると 正確な値とは言えないような 値が出てきました どこが間違っているのか説明していただければ 光栄です。 ちなみに データはこちら側にあります もっと詳しい情報が 提供できたら良いのですが 教えて Goo が 複数の写真を 送ることができないので このデータしかありません 申し訳ありません

質問者からの補足コメント

• ax^b で求まったaの値は5.977 またbの値は2.494でした。

    補足日時：2022/01/05 19:48

No.7 回答者： kamiyasiro 回答日時： 2022/01/06 16:48

#5です。

なお、両者の差は、添付図に示すように殆ど差がありません。強引に外挿すると、2.503乗の方が僅かに上を行きます。
差が小さい理由は、もともとの打点が誤差なくきれいに並んでいるからです。誤差が大きいと、どこをなるべく通そうとするかで曲線が違ってきます。

赤・・・a=6.015865, b=2.502520
青・・・a=5.984542, b=2.49476

No.6 回答者： kamiyasiro 回答日時： 2022/01/06 16:06

#5です。

訂正させて下さい。下記が正しいです。

#4では、
①変動係数一定仮定（これは#2さんの結果にほぼ一致）
②等分散仮定（これは#3さんの結果にほぼ一致）
③データ同化（a=1，b=1の初期値を与えて収束させた）→②に一致
を行っています。

ご質問者さんは①に近い結果を得られています。

No.5 回答者： kamiyasiro 回答日時： 2022/01/06 16:00

#4です。

答だけ投稿して、説明を怠り、申し訳ありません。

このような、非線形の関数フィッティングには、主に３つの方法があります。ご質問者は1)の方法を用いられたのであり、間違いではありません。

1)一般線形モデル・・・変数変換して、直線関係に持ち込んで、最小二乗法で解く。
2)一般化線形モデル・・・リンク関数のXの部分に線形式をブチ込んで、最尤法で解く。
3)データ同化・・・理論式のパラメータを変更しながら曲線をデータに近づける（代表的なものがレーベンバーグ・マルカート法）。

1)の欠点は、最小二乗法で行われる回帰だということです。それは、x重心y重心を通り、等分散の仮定で近似線を作ることと同義です。
そのため、ロジスティック曲線のような、０と１に漸近し、そこでは誤差は０である、というような曲線を引くことが出来ません。

2)はそれを解決しています。
今回の累乗関数（べき乗関数ではありません）は、０漸近していますので、(x,y)=(0,0)を通り、そこでは誤差は０であり、ｘが大きくなるにつれて誤差も増えると思われます。このようなモデルを、ユーザーが与えることができます。

#4では、
①等分散仮定（これは#3さんの結果にほぼ一致）
②変動係数一定仮定（これは#2さんの結果にほぼ一致）
③データ同化（a=1，b=1の初期値を与えて収束させた）→②に一致
を行っています。

ご質問者さんの求めた結果は、一般線形モデルで行って、たまたま②に近い結果が出たのであり、誤差モデルを意図して近似を行ったわけではありません。ここが理論的に弱い部分です。

ここを解決するのが一般化線形モデルなのです。

No.4 回答者： kamiyasiro 回答日時： 2022/01/06 10:35

Rのglm（一般化線形モデル）でやってみました。

もし、Rユーザーの方がいらっしゃったらご参考まで。

①このデータは０漸近するので、変動係数一定モデル（標準偏差がｘに比例して大きくなる）で回帰すると、(family = Gamma)
a=5.984542, b=2.49476

②０漸近とか無視して、等分散仮定で回帰すると、(family = gaussian)
a=6.015865, b=2.502520

③ちなみにＲのnlsで、原点を通る式で回帰すると、
a=6.01586, b=2.50252
②と③は等しくなります。

以下は実行したスクリプトと出力結果です。

> result <- glm(y ~ log(x), family = Gamma(link = log))
> summary(result)

Call:
glm(formula = y ~ log(x), family = Gamma(link = log))

Deviance Residuals:
Min 1Q Median 3Q Max
-0.075646 -0.005148 0.003347 0.010537 0.043496

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.78918 0.01542 116.0 3.41e-14 ***
log(x) 2.49476 0.01463 170.5 1.57e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.001035198)

Null deviance: 17.9195685 on 9 degrees of freedom
Residual deviance: 0.0085018 on 8 degrees of freedom
AIC: -40.069

Number of Fisher Scoring iterations: 3

> exp(result$coefficients[1])
(Intercept)
5.984542


> result <- glm(y ~ log(x), family = gaussian(link = log))
> summary(result)

Call:
glm(formula = y ~ log(x), family = gaussian(link = log))

Deviance Residuals:
Min 1Q Median 3Q Max
-0.031746 -0.010502 -0.002389 0.003920 0.034605

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.794400 0.002898 619.2 < 2e-16 ***
log(x) 2.502520 0.014168 176.6 1.18e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.0004324484)

Null deviance: 38.6226900 on 9 degrees of freedom
Residual deviance: 0.0034596 on 8 degrees of freedom
AIC: -45.313

Number of Fisher Scoring iterations: 3

> exp(result$coefficients[1])
(Intercept)
6.015865


> result <- nls(y ~ a * x^b + 0, data = X, start = list(a = 1, b = 1))
> summary(result)

Formula: y ~ a * x^b + 0

Parameters:
Estimate Std. Error t value Pr(>|t|)
a 6.01586 0.01743 345.1 < 2e-16 ***
b 2.50252 0.01417 176.6 1.18e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0208 on 8 degrees of freedom

Number of iterations to convergence: 6
Achieved convergence tolerance: 4.443e-08

No.3 回答者： stomachman 回答日時： 2022/01/05 22:46

ln(y[i]) = ln(a) + b ln(x[i]) + ε[i]

というモデルで
　　E = Σ(ε[i]^2)
を最小化したら
　　a = 5.987872013, b = 2.49578409, RMS = 0.022964813
となった。ここでRMSはroot mean square、すなわち √(E/10)のことです。
　よく合ってると思いますが、さらにもうちょっとだけ精度を上げられる。というのは本来なら
　　y[i] = a x[i]^b + δ[i]
というモデルで
　　F = Σ(δ[i]^2)
を最小化したいところなのに、対数を取ってfittingしたために、yが大きいほど誤差も大きい、という偏りが生じてしまっている。そこで、
　　y[i] = (a x[i]^b) (1 + δ[i]/(a x[i]^b))
の対数を取って
　　ln(y[i]) = ln(a) + bln(x[i]) + ln(1 + δ[i]/(a x[i]^b))
　　≒ ln(a) + bln(x[i]) + δ[i]/(a x[i]^b)
　　≒ ln(a) + bln(x[i]) + δ[i]/y[i]
だから、
　　y[i] ln(y[i]) ≒ ln(a) y[i] + b y[i] ln(x[i]) + δ[i]
である。そこで
　　y[i] ln(y[i]) = ln(a) y[i] + b y[i] ln(x[i]) + η[i]
というモデルで
　　G = Σ(η[i]^2)
を最小化するんです。つまりy[i]で重み付けしてfittingする。その結果は
　　a = 6.022866171, b = 2.505366152, RMS = 0.018365649

（以上の計算はExcelです。）

No.2 回答者： lon_donburi_dge 回答日時： 2022/01/05 20:42

> ax^b で求まったaの値は5.977 またbの値は2.494でした。

そんなに変な値？
代入すると、もとのyに結構合わない？
まあ、こちらで計算したら、a=5.988、b=2.496になりましたが。

計算過程をお教えいただけますか？

No.1 回答者： lon_donburi_dge 回答日時： 2022/01/05 16:33

> 正確な値とは言えないような値

a、bはいくつになったの？どこが間違っているのかと言われても
計算過程が分からないから、そこまではできないとしても、
a、bがいくつになったか次第では、何が間違った要因かわかるかも。