回答編集履歴

2

2018/10/04 11:05

投稿

tiitoi
tiitoi

スコア21956

test CHANGED
@@ -76,7 +76,7 @@
76
76
 
77
77
 
78
78
 
79
- ## 追記
79
+ ## 追記1
80
80
 
81
81
 
82
82
 
@@ -169,3 +169,91 @@
169
169
 
170
170
 
171
171
  ![イメージ説明](e8b555c0e450d809259f7cf52c3b1729.png)
172
+
173
+
174
+
175
+ ## 追記2
176
+
177
+
178
+
179
+
180
+
181
+ curve_fit() でも次のようにすることでパラメータを固定できました。
182
+
183
+ `stats.lognorm.fit()` のほうでは、fitting のアルゴリズムが異なるのか、データ数が少ない影響等により精度よく近似できませんでした。
184
+
185
+
186
+
187
+ ```python
188
+
189
+ import matplotlib.pyplot as plt
190
+
191
+ import numpy as np
192
+
193
+ from scipy import stats
194
+
195
+ from sklearn.metrics import mean_squared_error
196
+
197
+ from scipy.optimize import curve_fit
198
+
199
+
200
+
201
+ x = np.array([0, 4, 9, 14, 19, 24, 29, 34, 39, 44, 49, 54, 59,
202
+
203
+ 64, 69, 74, 79, 84, 89, 94, 99, 104, 109, 114, 119, 124,
204
+
205
+ 129, 134, 139, 144, 149, 154, 159, 164])
206
+
207
+
208
+
209
+ y = np.array([0., 0.17304493, 0.28618968, 0.50083195, 0.55407654,
210
+
211
+ 0.65058236, 0.73044925, 0.83527454, 0.87687188, 0.92845258,
212
+
213
+ 0.93510815, 0.95174709, 0.96006656, 0.9750416, 0.97670549,
214
+
215
+ 0.98169717, 0.98169717, 0.9843594, 0.98868552, 0.9906822,
216
+
217
+ 0.9906822, 0.99234609, 0.99234609, 0.99567388, 0.99567388,
218
+
219
+ 0.99567388, 0.99567388, 0.99567388, 0.99567388, 0.99567388,
220
+
221
+ 0.99733777, 1., 1., 1.])
222
+
223
+
224
+
225
+ # scipy.optimize.curve_fit を使うやり方
226
+
227
+ ######################################################
228
+
229
+ def cdf(x, a, b):
230
+
231
+ return stats.lognorm.cdf(x, s=a, loc=0, scale=b)
232
+
233
+ [s, scale], cov = curve_fit(cdf, x, y)
234
+
235
+ print('s={}, scale={}'.format(s, scale)) # s=0.8933602211719341, scale=14.750787612138023
236
+
237
+ # 近似した関数の結果
238
+
239
+ y_pred = cdf(x, s, scale)
240
+
241
+
242
+
243
+ # 描画する。
244
+
245
+ plt.plot(x, y, linestyle='--', marker='o', color='b', ms=2, label='data')
246
+
247
+ plt.plot(x, y_pred, linestyle='--', marker='o', color='g', ms=2, label='prediction')
248
+
249
+ plt.legend()
250
+
251
+ plt.show()
252
+
253
+
254
+
255
+ print(mean_squared_error(y, y_pred)) # 0.000730706490266028
256
+
257
+ print(y_pred[0]) # 0.0
258
+
259
+ ```

1

a

2018/10/04 11:05

投稿

tiitoi
tiitoi

スコア21956

test CHANGED
@@ -73,3 +73,99 @@
73
73
 
74
74
 
75
75
  ![イメージ説明](8d4ae870c9f93b672fab7f240017fa12.png)
76
+
77
+
78
+
79
+ ## 追記
80
+
81
+
82
+
83
+ ```
84
+
85
+ 対数正規分布の累積分布関数としているのに,xが0の時にyが0とならないのでどうにかしたいのです.
86
+
87
+ ```
88
+
89
+
90
+
91
+ 対数正規分布の定義域は 0 < x < ∞ なのに、その累積分布関数で cdf(0) = 0 とならないのはおかしいということですね。
92
+
93
+ cdf(0) = 0 とならない理由は scipy.stats.lognorm 関数に loc というシフトするパラメータを含んでいるからです。[リファレンス](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.lognorm.html) を参考にしてください。
94
+
95
+
96
+
97
+ なので、このまま fit() すると、loc も推定対象なので、近似したものでは cdf(0) = 0 とはなりません。scipy.stats.lognorm にも fit() という関数があるので、こちらを使うと、loc=0 と固定した状態で残りのパラメータを推定できます。
98
+
99
+
100
+
101
+ ```python
102
+
103
+ import seaborn as sns
104
+
105
+
106
+
107
+ from scipy import stats
108
+
109
+ from sklearn.metrics import mean_squared_error
110
+
111
+ from scipy.optimize import curve_fit
112
+
113
+
114
+
115
+ # 対数正規分布に従うサンプルを生成する。
116
+
117
+ shape, scale = 0.5, 1.
118
+
119
+ sample = stats.lognorm(s=0.5, loc=0, scale=1.).rvs(size=2000)
120
+
121
+ sns.distplot(sample, norm_hist=True, kde=False)
122
+
123
+ plt.show()
124
+
125
+
126
+
127
+ # loc は固定して、パラメータを推定する。
128
+
129
+ shape_pred, loc_pred, scale_pred = stats.lognorm.fit(sample, floc=0)
130
+
131
+ print('shape={}, loc={}, scale={}'.format(shape, loc, scale))
132
+
133
+ # shape=0.5094946936562328, loc=0, scale=1.008179589277641
134
+
135
+
136
+
137
+ x = np.linspace(0, 10, 100)
138
+
139
+ y = stats.lognorm.cdf(x, s=shape, loc=0, scale=scale)
140
+
141
+ y_pred = stats.lognorm.cdf(x, s=shape_pred, loc=loc_pred, scale=scale_pred)
142
+
143
+
144
+
145
+ # 描画する。
146
+
147
+ plt.plot(x, y, linestyle='--', marker='o', color='b', ms=2, label='true')
148
+
149
+ plt.plot(x, y_pred, linestyle='--', marker='o', color='g', ms=2, label='prediction')
150
+
151
+ plt.legend()
152
+
153
+ plt.show()
154
+
155
+
156
+
157
+ mse = mean_squared_error(y, y_pred)
158
+
159
+ print(mse) # 7.72357454555545e-06
160
+
161
+ print(y_pred[0]) # 0.00000000e+00 cdf(0) = 0 となっている。
162
+
163
+ ```
164
+
165
+
166
+
167
+ ![イメージ説明](6523c372a240c363576137d04c4a6cd8.png)
168
+
169
+
170
+
171
+ ![イメージ説明](e8b555c0e450d809259f7cf52c3b1729.png)