Rによる項目応答理論 p158 Newton-Raphson法を、PHPで書き直してみようという試みです。

NewtonRaphson.php

 
class NewtonRaphson {

  public array $a;
  public array $b;
  public array $u;
  public float $ε;

  public function __construct(array $a, array $b, array $u, float $ε) {

    $this->a = $a;
    $this->b = $b;
    $this->u = $u;
    $this->ε = $ε;

  }

  private function icc2PL(float $a, float $b, float $Θ) {

    return 1 / (1 + exp(- 1.7 * $a * ($Θ - $b)));;

  }

  private function LL(array $u, array $a, array $b, float $Θ) {

    $temp = [];

    for ($i = 0; $i < count($u); $i++) {

      $temp[] = ($u[$i] * log(self::icc2PL($a[$i], $b[$i], $Θ))) + ((1 - $u[$i]) * log(1 - self::icc2PL($a[$i], $b[$i], $Θ)));

    }

    return array_sum($temp);

  }

  private function dLL(array $u, array $a, array $b, float $Θ) {

    $temp = [];

    for ($i = 0; $i < count($u); $i++) {

      $temp[] = 1.7 * $a[$i] * ($u[$i] - self::icc2PL($a[$i], $b[$i], $Θ));

    }

    return array_sum($temp);

  }

  private function ddLL(array $a, array $b, float $Θ) {

    $temp = [];

    for ($i = 0; $i < count($a); $i++) {

      $temp[] = - pow(1.7, 2) * pow($a[$i], 2) * self::icc2PL($a[$i], $b[$i], $Θ) * (1 - self::icc2PL($a[$i], $b[$i], $Θ));

    }

    return array_sum($temp);

  }

  public function run() {

    // ある受験者の能力値 Θ 初期値を設定する
    // FROM Wainer et al 2007
    $Θ0 = log(array_sum($this->u) / (count($this->u) - array_sum($this->u)));

    echo '初期値 : ' . $Θ0 . "\n";

    $Θ1 = 0;        // 初期値
    $flag = false;  // フラグ

    while (!$flag) {

      $Θ1 = $Θ0 - (self::dLL($this->u, $this->a, $this->b, $Θ0) / self::ddLL($this->a, $this->b, $Θ0));

      if (abs($Θ1 - $Θ0) < $this->ε) {

        $flag = true;

      }

      $Θ0 = $Θ1;

    }

    echo '更新値 : ' . $Θ0 . "\n";

  }

}
 

index.php

 
require_once('NewtonRaphson.php');

$a = [1.0, 1.5, 0.7,  1.2]; // 識別力パラメタ
$b = [1.5, 0.5, 0.0, -1.0]; // 困難度パラメタ
$u = [1, 0, 1, 1];          // ある受験者の反応パタン
$ε = 0.0001;                // 収束基準用

$nr = new NewtonRaphson($a, $b, $u, $ε);
$nr->run();
 

出力結果

 
　初期値 : 1.0986122886681
　更新値 : 0.77643798426893