2016-08-01 53 views
2

我的模型基于以下增长函数:x_t + 1 = x_t + r * x_t *(1-x_t/K)增长率r由每年定义随机抽取(平均值= 0.2,标准偏差= 0.2)。我正在查看10000次运行70年后股票x的概率密度函数。当我计算概率密度函数下的面积时,对于0.1的标准偏差,其大致等于1,而对于0.2或0.3的标准偏差则不大。我用100个值为区域A创建了一个直方图。它在A = 2处显示非常高的峰值,但面积的值甚至达到8.为什么它不等于1?概率密度曲线下面积不等于1

runs=10000 

mean=0.2 
sd=0.2 

K=1 
x_0=0.4 
v=0.1 
t=70 

y=c() 
x=x_0 


for (j in 1:runs) { 

rand=rnorm(t,mean,sd) 

    for (i in 1:t) { 

    x=max(x+rand[i]*x*(x-v)*(1-x/K),0) 

    if(i==t) 
     y[j]=x 
    next} 

    x=x_0 
    next } 

    library(sfsmisc) 

    Dens=density(y) 

    f=approxfun(Dens$x, Dens$y) 

    h=c() 

    i=seq(0.9, 1, length.out=100000) 

     for (e in 1:length(i)) { 
     h[e]=f(i[e]) 
     next} 

    options(max.print=1000000) 

    h[is.na(h)]=0 

    area=sum(abs(h[-1]+h[-length(h)])/2*diff(i)) 

谢谢!

+1

你得到的是什么而不是1?如果它接近1 - 你期望什么?数值计算具有舍入误差。 –

+0

它很可能接近于2,但在8处的直方图中甚至有一些峰值,不能是舍入误差 –

+0

'library(sfsmisc)'的要点是什么? –

回答

0

您的数值积分存在很大的缺陷。

首先,这是密度的一个情节:

enter image description here

几乎所有的区域是一个非常狭窄的,令人难以置信的高穗下方。绝大多数的样本点都在这个峰值之外,因此在估计峰值面积方面做得不是很好。

如果评估print(Dens)你看到这一点:特别

Call: 
    density.default(x = y) 

Data: y (10000 obs.); Bandwidth 'bw' = 1.734e-06 

     x    y   
Min. :0.9922 Min. : 0.0 
1st Qu.:0.9942 1st Qu.: 0.0 
Median :0.9961 Median : 0.0 
Mean :0.9961 Mean : 439.1 
3rd Qu.:0.9981 3rd Qu.: 0.0 
Max. :1.0000 Max. :91131.2 

注意x值的所有都> 0.99,但您使用的是数值积分在区间[.9,1] 。因此,大部分数值积分都在密度范围之外。您正在使用(通过approxfun外推值在[0.9,0.99]的范围内,并且外推远远(相对来说)数据的范围。因此,90%的数值积分基于两个几乎为零的值(前两个y值)的不可靠外推。

您需要找到一些方法来重新调整数据。致电density本身是有问题的。它使用(作为默认值)512点,但如果你看看它们,只有几十个具有明显大于0的值。你试图通过大部分采样尾部而不是值来获得分布的图片发生重大事件的地方。