昨日ふと思いついて、プログラムを見直したら、バグが見つかりました。
それを直したら、ちゃんとポリコリック相関係数が算出できました！ばんざーい。
サイン関数を間違えていたのだ。サインといっても三角関数のsinではなくて、符号を取り出すsgn関数。数値計算の業界ではよく用いられるのだけどCの標準ライブラリには入っていない。簡単に作れるんだけど、面倒なのでネットで検索して引っかかったやつを使ったのだ。が、これがまずかった。返す型がint型だったので、double型の計算式の中に組み込むと0になっちゃうんだな。ご丁寧に式全体をキャストしてあったので、見つからなかったわけです。ちくしょー。
とりあえずうまく出たので、続いて方程式を解くアルゴリズムを考える。ニュートン・ラフソン法がいいのはわかっているんだけど、微分方程式を導出できないので、わかりやすいところで二分法にしようと。戸川先生の本を読むと、それより効率の良い割線法(Secant法)というのがあったので、それを使ってみる。瞬時にできたし、収束スピードも抜群。よしよし。
しかし、すぐ次の壁にぶつかるのであった。