teratail
回答率
85.48%
質問するログイン新規登録
トップに関する質問RのOptim関数の結果をニュートンラプソン法で同じ結果を算出したい

編集履歴

質問編集履歴

2

質問内容を修正した。

2019/11/02 22:56

投稿

51sep
51sep

スコア22

test CHANGED
File without changes
test CHANGED
@@ -54,7 +54,7 @@
54
54
 
55
55
  以下サイトのプログラムを拝見させて頂いていますが、なかなか確認できません。
56
56
 
57
- ニュートンラプソンループの前に、上記のようにmu=ax+bとし、
57
+ ニュートンラプソンループの前に、上記のようにlambda=exp(a*x+b)とし、
58
58
 
59
59
  a, bを算出したいのですが（この場合、0.05684214 2.43769624）、方法はありますでしょうか。
60
60
 

1

質問の例題を変更した。

2019/11/02 22:56

投稿

51sep
51sep

スコア22

test CHANGED
File without changes
test CHANGED
@@ -1,18 +1,18 @@
1
- RのOptim関数の結果をニュートンラプソン法で同じ結果を算出したい
1
+ Rのoptim関数の結果をニュートンラプソン法で同じ結果を算出したいと思います。
2
2
 
3
- Optim関数の結果をニュートンラプソン法で同じ結果を算出したいと思います。
3
+ 現状は以下のとおりです。ポワソン回帰です。
4
4
 
5
+ optim関数では以下プログラムで確認できています。
6
+
5
- 現状は以下のとおりです。どなたかご助言をお願いしてもよろしいでしょうか。
7
+ どなたかご助言をお願いしてもよろしいでしょうか。
6
8
 
7
9
 
8
10
 
9
- Optim関数では以下プログラムで確認できています。
10
-
11
11
  ```
12
12
 
13
- x<-c(1,1,1,1,1,2,2,2,2,2)
13
+ x<-c(1,2,10)
14
14
 
15
- y<-c(10,13,14,8,15,13,21,18,25,13)
15
+ y<-c(10,15,20)
16
16
 
17
17
  f=function(arg){
18
18
 
@@ -20,7 +20,7 @@
20
20
 
21
21
  b=arg[2]
22
22
 
23
- lambda=a*x+b
23
+ lambda=exp(a*x+b)
24
24
 
25
25
  -sum(y*log(lambda)-lambda-log(factorial(y)))
26
26
 
@@ -30,17 +30,17 @@
30
30
 
31
31
  #$`par`
32
32
 
33
- #[1] 6.002062 5.995965
33
+ #[1] 0.05684214 2.43769624
34
34
 
35
35
  #$value
36
36
 
37
- #[1] 27.00733
37
+ #[1] 7.150581
38
38
 
39
39
  #$counts
40
40
 
41
41
  #function gradient
42
42
 
43
- # 95 NA
43
+ # 57 NA
44
44
 
45
45
  #$convergence
46
46
 
@@ -52,98 +52,14 @@
52
52
 
53
53
  ```
54
54
 
55
- 以下サイトのプログラムを使用させて頂いています。
55
+ 以下サイトのプログラムを拝見させて頂いていますが、なかなか確認できません。
56
56
 
57
- このニュートンラプソンループの前の“poisson.lik”に、上記のようにxの情報を加えてmu=ax+bとし、
57
+ ニュートンラプソンループの前に、上記のようにmu=ax+bとし、
58
58
 
59
- 6.002062、5.995965、27.00733を算出したいと思います。
59
+ a, bを算出したいのですが（この場合、0.05684214 2.43769624）、方法はありますでしょうか。
60
60
 
61
- 単純に上記の“f”のように中身を変えても、うまく実行できません。
62
-
63
- どなたかご助言をお願いしてもよろしいでしょうか。
61
+ 具体的な参考例をご存知の方、ご助言をお願いしてもよろしいでしょうか。
64
62
 
65
63
 
66
64
 
67
65
  [Programming Newton Raphson in R for Maximum Likelihood estimation](https://stackoverflow.com/questions/42683458/programming-newton-raphson-in-r-for-maximum-likelihood-estimation)
68
-
69
- ```# NEWTON-RAPHSON METHOD
70
-
71
- #generate the data
72
-
73
- lambda <- 3.2
74
-
75
- ydata <- rpois(500, lambda)
76
-
77
-
78
-
79
- #declare the log likehood function
80
-
81
- poisson.lik <- function(mu = ""){
82
-
83
- n <- length(ydata)
84
-
85
- loglik <- -n * mu - sum(log(factorial(ydata))) + log(abs(mu)) * sum(ydata)
86
-
87
- return(-loglik)
88
-
89
- }
90
-
91
-
92
-
93
- #creating the Newton Raphson loop
94
-
95
- NR <- function ( mu = "", initval = "", f = "", stoptol = 1e-05, imax = "") {
96
-
97
- i = 0
98
-
99
- h = .1
100
-
101
- mu0 = initval - 0.1
102
-
103
- mu1 = initval
104
-
105
- df.dx <- double(1) #predeclare your variable types
106
-
107
- while ( abs(f(mu) - f(mu1)) > stoptol && i <= imax ) {
108
-
109
- mu0 <- mu1
110
-
111
- df.dx <- (f(mu0 + h) - f(mu0)) / h
112
-
113
- mu1 <- mu0 + (mu - mu1) / abs(df.dx)
114
-
115
- i <- i + 1
116
-
117
- }
118
-
119
- return(list("nstep" = i, "Final" = mu1, "fctval" = f(mu1)))
120
-
121
- }
122
-
123
-
124
-
125
- mu <- 3.2
126
-
127
- initval <- 1
128
-
129
- f <- poisson.lik
130
-
131
- stoptol <- 0.0001
132
-
133
- imax <- 10^4
134
-
135
- NR(mu, initval, f, stoptol, imax)
136
-
137
- #$`nstep`
138
-
139
- #[1] 623
140
-
141
- #$Final
142
-
143
- #[1] 3.20002
144
-
145
- #$fctval
146
-
147
- #[1] 991.2795
148
-
149
- ```