Hubery-Lee

蒹葭苍苍，白露为霜

0%

数据拟合

发表于 分类于 Valine：
介绍matlab,root,python如何使用进行自定义函数进行数据拟合

数据拟合

[toc]

1.问题描述

手里有三组数据，用来分析探测器散射体材料厚度与探测器效率之间的关系。需要对三组数据进行画散点图，并进行拟合。拟合函数采用如下形式
$$
\eta = 1- exp(-\mu x)
$$
准备采用三种方法进行数据拟合，数据分析软件root、matlab和python.

2.三种数据分析软件的数据拟合实现

2.1 matlab的数据拟合实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
%% 读入数据x,y
x=5:5:100;
y1=textread('dataout1.txt'); //textread读入数据
y2=textread('dataout2.txt');
y3=textread('dataout3.txt');

%% 自定义函数形式数据拟合
myfittype = fittype('1. - exp(a*x)',...
    'dependent',{'y'},'independent',{'x'},...
    'coefficients',{'a'});

myfit1 = fit(x',y1,myfittype);
plot(myfit1,'r-',x,y1,'r*');
hold on
myfit2 = fit(x',y2,myfittype);
plot(myfit2,'b-',x,y2,'bo');
hold on
myfit3 = fit(x',y3,myfittype);
plot(myfit3,'g-',x,y3,'gx');

fittedData

2.2 root中的数据拟合实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
//参考 root tutorial/graphs/zdemo.C  /fit/FittingDemo.C

#include "iostream"
#include "fstream"
#include "vector"
#include "cstring"
#include "string.h"

using namespace std;

const int NMAX=20;
double xx[NMAX],yy[NMAX];

class ReadData
{
public:
 
  ReadData(string fname)
  {
    filename=fname;
  };
  ~ReadData(){};

  void GetData();
private:
  string filename;
};

void ReadData::GetData()
{
  ifstream input(filename);
  double a,b;
  vector<double> x,y;
  int i=0;
  while(!input.eof())
    {
      i++;
      a=i*5;
      input>>b;
      x.push_back(a);
      y.push_back(1.0-(b/1.e6));
    }
  input.close();
  i=i-1;

  for(int j=0;j<i&&j<NMAX;j++)
    {
      xx[j]=x[j],yy[j]=y[j];
    }
}



void plot()
{

  ReadData a1("./EJ200-100/data.txt");
  a1.GetData();

  TGraph *g1 =new TGraph(NMAX,xx,yy);
  g1->SetLineColor(38);
  g1->SetMarkerColor(kBlue);
  g1->SetMarkerStyle(21);
  g1->SetMarkerSize(1.1);
  g1->Draw("AP");
  
  ReadData a2("./EJ200-140/data.txt");
  a2.GetData();

  TGraph *gr2 = new TGraph(NMAX,xx,yy);
  gr2->SetLineColor(38);
  gr2->SetMarkerColor(kRed);
  gr2->SetMarkerStyle(29);
  gr2->SetMarkerSize(1.5);
  gr2->Draw("P");

  ReadData a3("./EJ200-180/data.txt");
  a3.GetData();

  TGraph *gr3 = new TGraph(NMAX,xx,yy);

  gr3->SetLineColor(38);
  gr3->SetMarkerColor(6);
  gr3->SetMarkerStyle(8);
  gr3->SetMarkerSize(1.1);
  gr3->Draw("P");

  //////////////////
  ////
  //      elta = 1-exp(a*x) //拟合函数形式

  TF1 *fun =new TF1("#elta = 1-exp(A*x)","1.-exp([0]*x)",0,100);
  fun->SetLineColor(kBlue);fun->SetLineStyle(2);
  g1->Fit(fun);

  fun->SetLineColor(kRed);fun->SetLineStyle(2);
  gr2->Fit(fun);

  fun->SetLineColor(6);fun->SetLineStyle(2);
  gr3->Fit(fun);

  //////////////////
  //
  //  legend
  TLegend *lg=new TLegend(0.6,0.65,0.88,0.85);
  lg->SetTextFont(62);
  lg->SetTextAlign(12);
  lg->SetTextSize(0.04);
  lg->AddEntry(g1,"100keV","p");
  lg->AddEntry(gr2,"140keV","p");
  lg->AddEntry(gr3,"180keV","p");
  lg->AddEntry(fun," #eta = 1-exp(#mu *x)","l");
  lg->Draw();
}

c1

2.3 python中的数据拟合实现

python 中做数据分析常用的几个包matplotlib,numpy,scipy,pandas

  • scipy中的curve_fit函数
  • lmfit中的fit函数

下面从这两种包来拟合这两组数据。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
# -*- coding: utf-8 -*-
# @Time    : 2020/6/6 17:52
# @Author  : Hubery-Lee  
# @Email   : hrbeulh@126.com

##   case 1: curve_fit
##   case 2: lmfit
##
# Header
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from lmfit import Model


# Define function
def func(x, a):
    return 1. - np.exp(a * x)


# Read data from text
x = np.linspace(5, 100, 20)
y1 = np.loadtxt('dataout1.txt')
y2 = np.loadtxt('dataout2.txt')
y3 = np.loadtxt('dataout3.txt')
#  plot data
plt.plot(x, y1, 'bo', label='100 keV')
plt.plot(x, y2, 'r*', label='140 keV')
plt.plot(x, y3, 'gx', label='180 keV')

## case 1:
## https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html
# Fit for the parameters a of the function func:
popt1, pcov1 = curve_fit(func, x, y1)
# popt  # output: array([ 2.55423706, 1.35190947, 0.47450618])
plt.plot(x, func(x, *popt1), 'b-',
         label='100 keV fit: a=%5.3f' % tuple(popt1))

popt2, pcov2 = curve_fit(func, x, y2)
# popt  # output: array([ 2.55423706, 1.35190947, 0.47450618])
plt.plot(x, func(x, *popt2), 'r-',
         label='140 keV fit: a=%5.3f' % tuple(popt2))

popt3, pcov3 = curve_fit(func, x, y3)
# popt  # output: array([ 2.55423706, 1.35190947, 0.47450618])
plt.plot(x, func(x, *popt3), 'g-',
         label='180 keV fit: a=%5.3f' % tuple(popt3))

# # In the case of parameters a need be constrainted
# # Constrain the optimization to the region of
# # 0 <= a <= 3, 0 <= b <= 1 and 0 <= c <= 0.5
# popt, pcov = curve_fit(func, xdata, ydata, bounds=(0, [3., 1., 0.5]))
# popt  # output: array([ 2.43708906, 1. , 0.35015434])
# plt.plot(xdata, func(xdata, *popt), 'g--',
#          label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))

## case 2:
## https://lmfit.github.io/lmfit-py/model.html
# Fitting
gmodel = Model(func)
result = gmodel.fit(y1, x=x, a=-0.02)  # Fit from initial values (5,5,1)
print(result.fit_report())

# Plot
#plt.plot(x, y1, 'bo', label='raw data')
#plt.plot(x, result.init_fit, 'b--', label='init_fit')
#plt.plot(x, result.best_fit, 'k--', label='best_fit')
plt.plot(x, result.best_fit, 'k--', label='100 keV lmfit')


# draw option
# Labels
plt.title(r'Fitting Function $\eta = 1.0-exp(a*x)$')
plt.xlabel('x/mm')
plt.ylabel(r'$\eta$')  # 公式的添加 latex风格 https://matplotlib.org/tutorials/text/mathtext.html
plt.legend()
leg = plt.legend()  # remove the frame of Legend, personal choice
leg.get_frame().set_linewidth(0.0)  # remove the frame of Legend, personal choice
# leg.get_frame().set_edgecolor('b')  # change the color of Legend frame
# plt.show()

# Export figurey
# plt.savefig('fit1.eps', format='eps', dpi=1000)
# plt.savefig('fit1.pdf', format='pdf', dpi=1000, figsize=(8, 6), facecolor='w', edgecolor='k')
plt.savefig('myfit.jpg', format='jpg', dpi=1000, figsize=(8, 6), facecolor='w', edgecolor='k')

myfit

Donate comment here.

欢迎关注我的其它发布渠道