如何在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]