直方图中的概率密度函数(gnuplot)

问题描述 投票:0回答:1

我一直在 gnuplot 中绘制直方图,但我还想在直方图本身之上包含概率密度分布。目前我可以生成下图,但概率密度函数本身的尺度似乎完全偏离了。

来自“DZ_ref-CCSD_error”的数据 - https://pastebin.com/xX1MUd8i

附代码:

reset
stats 'DZ_ref-CCSD_error' using 1
 n=STATS_records #number of intervals
 max=STATS_max #max value
 min=STATS_min #min value
 width=(max-min)/n #interval width
 binwidth=6
#+width/2.0 to center graph
 bin(x,width)=width*floor(x/width) + width/2.0
 set term png
 set output "histogram_out.png"
 set xtics -100,50,100
 set boxwidth binwidth
 set mxtics 5
 set style fill solid 0.5 
 set tics out nomirror
 set tics font ", 12"
 set ytics 0,0.25,0.5
 set mytics 5
 set yrange [0:0.50]
 set xrange [-100:100]
 gauss(x)=a/(sqrt(2*pi)*sigma)*exp(-(x-mean)**2/(2*sigma**2))
 fit gauss(x) 'DZ_ref-CCSD_error' via a, sigma, mean
 plot 'DZ_ref-CCSD_error' u (bin($1, binwidth)-0.5*binwidth):(1./n*binwidth) smooth fnormal w boxes lc           rgb     "green" notitle, \
  gauss(x) w lines ls 2 lw 2

概率密度不正确的直方图

如果有人能指出我的错误在哪里或对这个问题添加任何见解,我将不胜感激!

gnuplot histogram probability-density
1个回答
0
投票

一个简单的错误是您试图拟合原始数据,但实际上需要拟合分箱的直方图数据。 因此,简单的解决方案:将直方图数据写入数据块表,例如

$Histo
并对这些数据进行高斯拟合。

检查以下完整(稍微重写)的示例:

数据:

SO79153881.dat

-26.426
-0.062
-2.207
-0.302
-2.128
-5.268
7.564
-5.959
0.016
1.069
5.012
-2.561
16.145
27.64
-0.049
-6.315
0.0597
0.047
122.326
0.354
0.003
122
3.016
1.122
0.001
0.172
0.906
-13.012
9.815
-9.193
-0.041
0.001
-10.143
-4.746
17.549
-26.467
-0.015
0.254
6.102
-3.09
12.202
-9.759
1.793
25.296
-2.096
1.751
0.35
2.916
-0.43
0.003
0.395
3.813
29.021
3.498
-5.013
0.401
0.89
1.543
-1.611
1.159
0.055
2.949
6.737
-5.935
0
-0.015
-7.406
0.471
1.45
2.37
-4.509
1.473
0.002
0.941
0.285
-0.029
1.607
-0.729
-20.279
19.272
2.189
0.0048
4.144
0.181
0.083
0.162
-2.71
1.699
0.084
0.524
0.19
-0.011
0.08
3.562
26.203
7.459
-6.999
1.833
2.926
3.28
2.247
30.148
-0.098
26.572
-15.62
-0.119
0.248
-0.235
-13.173
11.535
8.358
-0.855
2.007
-2.621
7.28
0.808
-0.484
28.104
-0.021
2.023
10.761
1.195
-1.966
9.326
0.006
1.137
0.268
0.303
-4.066
12.69
-0.013
0.019
0.579
1.675
-0.076
0.006
-8.134
-1.428
-0.854
0.408
-0.06
7.55
0.215
0.647
0.001
-0.358
-12.749
-0.009
-35.656
-4.529
-11.391
0.008
-0.012
-0.132
3.561
-22.863
6.031
4.271
-0.12
6.817
10.17
-0.089
61.425
1.538
0.045
-0.184
-15.299
0.195
0.039
-0.087
-6.891
0.002
0.475
-0.118
-0.977
-10.807
-0.066
-2.457
0.02
-0.002
4.929
-1.562
-0.016
0.826
-0.047
1.202
0.073
1.383
-0.245
-5.766
-0.002
1.34
0.223
-1.505
-0.737
-0.013
-0.012
-0.018
-0.006
-14.78
0.368
-0.434
3.16
8.424
0.132
1.64
-0.22
35.201
0.458
0.005
0.847
27.457
-0.032
0.029
0.216
11.494
-0.035
0.026
0.767
-3.385
-0.871
1.97
1.222

脚本:

### fit gaussian to histogram
reset session

FILE = "SO79153881.dat"

stats FILE u 1
n        = STATS_records     # number of intervals
max      = STATS_max         # max value
min      = STATS_min         # min value
width    = (max-min)/n       # interval width
binwidth = 6
bin(x,width) = width*floor(x/width) + width/2.0

set xtics -100,50,100
set boxwidth binwidth
set mxtics 5
set style fill solid 0.5 
set tics out nomirror
set tics font ", 12"
set ytics 0,0.25,0.5
set mytics 5
set yrange [0:0.50]
set xrange [-100:100]
set samples 500

set table $Histo
     plot FILE u (bin($1, binwidth)-0.5*binwidth):(1./n*binwidth) smooth fnormal
unset table

gauss(x) = a/(sqrt(2*pi)*sigma)*exp(-(x-mean)**2/(2*sigma**2))
fit gauss(x) $Histo u 1:2 via a, sigma, mean

plot $Histo u 1:2 w boxes lc "green" notitle, \
     gauss(x) w lines ls 2 lw 2
### end of script

结果:

enter image description here

© www.soinside.com 2019 - 2024. All rights reserved.