scipy.optimize.approx_fprime

先日、フィッシャー情報量を対数尤度の2次の導関数から計算してみようと思い、ChatGPTに「pythonで数値微分するコードの例を出して」と尋ねてみました。

するとscipyの approx_fprime という関数が提案され、2次の導関数については approx_fprime を2回再帰的に適用する方法が提案されました。

これは例えば以下のように書くことができます。

import numpy as np
from scipy.optimize import approx_fprime
EPSILON = np.sqrt(np.finfo(float).eps)

def f(x):
    """微分したい対象の関数"""
    return np.sin(x)

def df(x):
    """1次の導関数"""
    return approx_fprime(x, f=f, epsilon=EPSILON)[0]

def ddf(x):
    """2次の導関数"""
    return approx_fprime(x, f=df, epsilon=EPSILON)[0]

ただ、これplotしてみると2階微分のほうはギザギザした曲線になるみたいです。

どうしてこうなるのかは approx_fprime の中身を知らないのでよくわからないのですが、とりあえず別の方法を探すことにしました。

numdifftools

numdifftools は数値微分をするためのパッケージで、 numdifftools.Derivative([対象の関数], n=[導関数の次数]) といった書き方で簡単に導関数を生成できます。

import numpy as np
import numdifftools as nd

def f(x):
    """微分したい対象の関数"""
    return np.sin(x)

df = nd.Derivative(f, n=1)
ddf = nd.Derivative(f, n=2)

これでプロットすると