Inverf.java

/*
 * $Id: Inverf.java,v 1.17 2008/03/03 15:18:35 koga Exp $
 * 
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.matrix;

import org.mklab.nfc.matrix.DoubleMatrix;


/**
 * 逆誤差関数を表すクラスです。
 * 
 * <p> Inverse Error function
 * 
 * @author koga
 * @version $Revision: 1.17 $
 * @see org.mklab.tool.matrix.Erf
 */
public class Inverf {

  /**
   * <code>x</code>の全ての成分の逆誤差関数を計算します。 精度は、<code>1E-5</code>です。
   * 
   * <p>誤差関数は、正規分布の確率分布関数を <code>0</code>から<code>x</code>まで積分する
   * 
   * <pre><code> erfnc(x) = 2/sqrt(PI) integral(0,x) exp(-t&circ;2) </code></pre>
   * 
   * です。
   * 
   * @param x データ
   * @return 逆誤差関数 (inverse error function)
   */
  public static DoubleMatrix inverf(DoubleMatrix x) {
    double tolerance = 0.00001;
    return inverf(x, tolerance);
  }

  /**
   * <code>x</code>の全ての成分の逆誤差関数をけいさんします。 ただし、精度は<code>tol</code>です。
   * 
   * @param x データ
   * @param tolerance 許容誤差
   * @return 逆誤差関数 (inverse error function)
   */
  public static DoubleMatrix inverf(DoubleMatrix x, double tolerance) {
    double t0 = 1 / Math.sqrt(2 * Math.PI);
    DoubleMatrix b0 = x.createOnes(x.getRowSize(), x.getColumnSize());
    DoubleMatrix b = x.createZero(x.getRowSize(), x.getColumnSize());

    // Newton-Raphson iteration
    while (b.subtract(b0).absElementWise().compareElementWise(".>", tolerance).anyTrue()) { //$NON-NLS-1$
      b0 = b;
      DoubleMatrix f = Erf.erf(b0).subtract(x).divide(2);
      DoubleMatrix fd = b0.multiplyElementWise(b0).multiply(0.5).unaryMinus().expElementWise().multiply(t0);
      b = b0.subtract(f.divideElementWise(fd));
    }
    return b;
  }

}