2016-09-28 122 views
-1

在R,I需要计算条件期望F(z)的= E [X | X < Z],其中X分布根据参数分布(即,对数正态分布)。计算期望(比如对数正态)中的R

为了计算例如F(2)我已经做了以下内容:

zz <- rlnorm(1000,meanlog=.7,sdlog=.5) 
mean(zz[zz<2]) 

不过,我不知道是否有一个更直接的方式,不需要生成样本。

回答

1

您正在查看truncated distribution。将x * f(x)积分为(-Inf, z),然后将该积分除以F(z)。 [f(x)是无条件的PDF; F(x)是无条件CDF。]

## integrand 
f <- function(x, mu, sigma) x * dlnorm(x, mu, sigma) 

## conditional expectation 
g <- function(z, mu, sigma) { 
    int <- integrate(f, lower = -Inf, upper = z, mu = mu, sigma = sigma) 
    int$value/plnorm(z, mu, sigma) 
    } 

## theoretical value 
g(2, 0.7, 0.5) 
# [1] 1.401472 

## sample estimate 
set.seed(0) 
zz <- rlnorm(1000,meanlog=.7,sdlog=.5) 
mean(zz[zz<2]) 
# [1] 1.40316 

我已经刨去写乳胶一行或两行说明了为什么我们需要一个整体如上,但它看起来像维基百科的链接是足够的信息。


For some reason, I am not able to plot the resulting function g . plot(1, g(1:10,0.7,0.5)) is giving an error.

要绘制,你需要使它成为一个量化的功能第一g。有一些关于绘制积分的帖子,如R plotting integral。下面是我们可以做的:

vg <- Vectorize(g, vectorize.args = "z") 
plot(1:10, vg(1:10, 0.7, 0.5), type = "l") 

enter image description here

+0

这似乎是一个很好的解决方案。由于某种原因,悬停,我无法绘制结果函数g。如果我做'z = 1:10',然后'plot(z,g(z,2,3))',那么结果图是不正确的。事实上,如果我重新定义'g',以便它不被'plnorm(...)分割',我得到以下错误:'xy.coords(x,y,xlabel,ylabel,log)中的错误:'x'和'y'长度不同' – Massimo2013

0

宋哲元通过的回答启发,做了一点点研究上的功能,其中的条件概率密度为截断PDF的条件期望。

据我,mean(zz[zz < a])在有条件的宇宙X <一个条件期望,因为这是用于生成ZZ值是原来的对数正态分布的PDF而不是的PDF有条件截短的pdf

为了计算条件期望我们必须使用截断PDF截断分布和不是原来的对数正态分布得出样本,然后计算期望。

如可以看到的,的mean(zz[zz < a])值总是从条件期望不同使用期望具有条件(截断)PDF计算,差值随着一个增加(任何直觉为什么?)。

# compute the truncated pdf with x < a 
tr.pdf <- function(x, a, m, s) (x < a) * (dlnorm(x, m, s)/plnorm(a, m, s)) 

expect.f <- function(x, a, m, s) x * tr.pdf(x, a, m, s) 

cond.expect.f <- function(a, m, s) { 
    return(integrate(expect.f, lower = -Inf, upper = a, a = a, m = m, s = s)$value) 
} 

m <- .7 
s <- .5 
curve(tr.pdf(x, a=2, m, s), 0, 5, col='red', ylab='y') 
curve(tr.pdf(x, a=2.5, m, s), 0, 5, col='green', add=TRUE) 
curve(tr.pdf(x, a=3, m, s), 0, 5, col='blue', add=TRUE) 
curve(dlnorm(x, m, s), 0, 5, add=TRUE) 

enter image description here

n <- 100000 
zz <- rlnorm(n,meanlog=m,sdlog=s) 

a <- 2 
mean(zz[zz<a]) 
#[1] 1.404279 
cond.expect.f(a, m, s) 
#[1] 1.401472 

a <- 2.5 
mean(zz[zz<a]) 
#[1] 1.622174 
cond.expect.f(a, m, s) 
#[1] 1.617784 

a <- 3 
mean(zz[zz<a]) 
#[1] 1.794217 
cond.expect.f(a, m, s) 
#[1] 1.787772 

对这个有什么想法?

+0

在我看来,你所绘制的密度函数在相关区间内都是相同的(除了它们被归一化以便它们之间的区域为1)。您在值中看到的差异只是“精确”值与样本值之间的差异(如果您生成新样本,则可以看到样本均值可以高于或低于计算值) – Massimo2013

+0

@Massimo:I已经使用了相当高的值,这就是为什么平均值的波动将会更小。但是我的观点是不同的:当你计算平均值(zz [zz

+0

正如我写的,在相关区间(z Massimo2013