如何在java中擬合正態分布

前言

最近在工作中需要擬合高斯曲線，在python中可以使用 scipy,相關代碼如下：

#!/usr/bin/env python 
# -*- coding=utf-8 -*-
%matplotlib inline

import numpy as np
import pylab as plt


from scipy.optimize import curve_fit

x = range(10)
y = [25,68,144,220,335,199,52,14,5,2]

def gaussian2(x,*param):
    return param[0]*np.exp(-np.power(x - param[1], 2.) / (2 * np.power(param[2], 2.)))


def use_fit_gaussian2():
    try:
        popt,pcov = curve_fit(gaussian2,x,y,p0=[1,1,1])
    except Exception,e:
        print e
        return

    print popt
    
    yhat = gaussian2(x, *popt)                         # or [p(z) for z in x]
    ybar = np.sum(y)/len(y)          # or sum(y)/len(y)
    ssreg = np.sum((yhat-ybar)**2)   # or sum([ (yihat - ybar)**2 for yihat in yhat])
    sstot = np.sum((y - ybar)**2)    # or sum([ (yi - ybar)**2 for yi in y])
    error = ssreg / sstot
    print ssreg, sstot, error

    plt.plot(x,y,'b+:',label='data', linewidth =3)
    plt.plot(x,gaussian2(x,*popt),'ro:',label='fit', linewidth =3)
    plt.legend()
    plt.show()

use_fit_gaussian2()

生成的結果如下圖所示：

java ?

由於線上用的java，所以需要使用java實現，需要使用到 apache 的 commons-math3 jar包

 <dependency>
     <groupId>org.apache.commons</groupId>
     <artifactId>commons-math3</artifactId>
     <version>3.6.1</version>      
 </dependency>

• 代碼

import org.apache.commons.math3.fitting.GaussianCurveFitter;
import org.apache.commons.math3.fitting.WeightedObservedPoints;

/**
 * Created by xingxing.dxx on 2016/11/14.
 */
public class CurveFittingTest {

    public static void main(String[] args) {
        WeightedObservedPoints obs = new WeightedObservedPoints();
        obs.add(0, 25);
        obs.add(1, 68);
        obs.add(2, 144);
        obs.add(3, 220);
        obs.add(4, 335);
        obs.add(5, 199);
        obs.add(6, 52);
        obs.add(7, 14);
        obs.add(8, 5);
        obs.add(9, 2);
        
        double[] parameters = GaussianCurveFitter.create().fit(obs.toList());
        for (double i : parameters) {
            System.out.println(i);
        }
    }
}

最開始測試的時候非常完美，可是馬上就悲劇了，在我運行如下case的時候，拋了一個錯

[0, 0, 0, 0, 4, 1, 0, 2, 0]

Exception in thread "main" org.apache.commons.math3.exception.ConvergenceException: illegal state: unable to perform Q.R decomposition on the 9x3 jacobian matrix

瞬間懵逼了，有沒有。然后開始找錯誤是哪報的，發現了報錯的代碼,

 if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
           throw new ConvergenceException(LocalizedFormats.UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN,
                                                   nR, nC);
 }

看着代碼貌似是矩陣求逆的時候，矩陣中有1.0/0.0的情況，猜測是不是輸入的數據都大於0就ok了，果然，下面的測試樣本就ok了

[0.01, 0.01, 0.01, 0.01, 4.01, 1.01, 0.01, 2.01, 0.01]