Scipyのoptimize.fmin関数を使った階段法における心理物理曲線

最近，初めて階段法を使った心理物理実験をやりました．
実験の後半は僕も被験者も意識が朦朧としてました...笑
ただ，プロトコルしっかり練って，それに対応したプログラム組めば実験者側はかなり楽になると思います．
被験者側はしっかり休憩挟む以外，楽する方法は無さそうですね笑
心理物理実験の後はその結果に対してカーブフィッティングを行うことで，PSE（主観的等価点），75％閾値などを算出して基準に対する比較刺激のバイアスがあるかどうかを調べていきます．
カーブフィッティングの方法については下記の梶本先生のページが一番わかりやすいと思います．
恒常法のS字カーブフィッティングについて
せっかくMatlabで書いて頂いているのですが，それを参考にして今回はPythonで書いてみました．

まぁこれといってコード自体に大きな違いはないのですが，標準正規累積分布関数のμとσを探索する際，Matlabのfminserachの代わりにScipyのoptimize.fmin関数を使って探索を行います．

Scipy.optimize.fmin()

リファレンスはこちらになります．
scipy.optimize.fmin — SciPy v1.0.0 Reference Guide
簡単に今回使った中身の変数を紹介していくと，

fmin(最小値探索を行う目的関数, 最適化を行う変数, args=目的関数に入力する定数)

という感じになります．
梶本先生のページのコード内における変数を使って書くと，

 fmin(err, numpy.array([1,cs.mean()]), args = (cs,p)) 

こうなります（Pythonの記述法に沿ったものに少し書き換えてます）．
errは

import math
import scipy as sp

def err(x,cs,p):
    sigma,mu = x
    a = p - 0.5 - 0.5 * sp.erf((cs-(mu))/(math.sqrt(2)*sigma))
    return sum((a**2))

で定義されます．

ただ，今回は階段法なので比較刺激によっては繰り返した回数が異なるため，比較刺激におけるデータの質がそれぞれ異なります．
そこで，変数を探索するための目的関数をただの最小二乗法ではなく，重み付き最小二乗法にします．
jp.mathworks.com
（Matlabのリファレンスですが，これが一番わかりやすかったです）
重み付きにすることで，カーブフィッティングをする際に質の高いデータ（閾値付近）が優先的に考慮されるようになります．
重み付きにしたからといって特に記述が複雑になるわけではなく，fminにおいては

 fmin(err, numpy.array([1,cs.mean()]), args = (cs,p,w)) 

argsに重み配列であるwを新たに追加し，errにおいては

import math
import scipy as sp

def err(x,cs,p,w):
    sigma,mu = x
    a = p - 0.5 - 0.5 * sp.erf((cs-(mu))/(math.sqrt(2)*sigma))
    return sum(w*(a**2))

返り値に重みを乗じるだけで済みます．

結果はこのようになります．

青線が普通の最小二乗法を使ったカーブフィッティング，赤線が重み付き最小二乗法を使ったカーブフィッティングになっています．
比較刺激1.5あたりで外れ値があり，青線はその外れ値に少し引っ張られていますが赤線はしっかりと閾値付近（繰り返し数が多い部分）で立ち上がっています．
これはあくまで例で，実験中に外れ値はめったに計測されないのですが，階段法だと上記したように同じ刺激を何度も繰り返すことになるので
外れ値対策というよりはしっかりと人間の弁別閾を見極められるのが重み付き最小二乗法のメリットだと思います．
コード全体についてはGithubの方へアップしておりますのでそちらを参考にしてください．
github.com

今回記述しませんでしたが，コード内では収束したかどうかを調べるためのグラフも出力しております．
そちらに関しては，下記のページを参考にして書いております．
qiita.com

間違い等あれば是非指摘してください．