0

我在使用 MATLAB 时遇到以下问题：

令 Z 服从对数正态分布，使得 ln Z 具有均值 m 和方差 w。设 eta 为负数，ca 为正常数。

我正在尝试计算期望值（让 I(Z<=c) 表示集合的指示函数 (Z<=c)）

E[Z^(eta+1) I(Z<=c)] = (1/sqrt(w)) 积分_0^cx^(eta) phi((ln x - m)/sqrt(w)) dx,

其中 phi() 表示标准正态随机变量的概率分布函数。

我做的第一件事是模拟 10.000 次 Z 试验，将值 >c 的向量的条目设置为 0，提高到 (eta+1) 的幂，然后计算平均值。这应该给我预期值的 MC 估计。

ST = random('Lognormal', m,w_sq,10000,1);
hlp = zeros(10000,1);
hlp(ST<=2) = ST(ST<=2);
hlp(hlp>0) = hlp(hlp>0).^(eta1+1); % 0^(eta1+1) gives infinity
mean(hlp)

对于积分，我使用了以下代码

tmpp = integral(@(x) x.^(eta1) .* normpdf((log(x)-m)/sqrt(w_sq),0,c);
tmpp / sqrt(w_sq(1))

不幸的是，这些程序导致了完全不同的结果，尽管在数学上它们应该是相同的。

整个事情是更大代码的一部分，对我来说使用积分版本会更方便。最初我尝试使用 MC 模拟仔细检查，然后发现一定有问题......

有人可以帮忙吗？

4

1 回答 1

0

对于第一段代码，

我猜你有意外结果的原因是当你在 上进行计算时hlp，你试图避免 0 值，因为这0^eta 会爆炸 - 这不是想要的结果，所以你只需放弃它。但在最后一步，mean(hlp)将获取 arrayhlp中的所有值，包括那些 0。尝试这个：

mean(hlp(hlp>0)).

我的结果在 10,000 点模拟时大约下降 2.x，在 1,000,000 点模拟时下降 2.3x ~ 2.4x。

我错了。您的问题是您使用的积分太少。尝试 10,000,000 积分，您会满意 :)

其次，我在理解您定义变量的方式时遇到了问题。（我没有足够的声誉来添加评论，所以我把它放在这里。）中的“sq”是w_sq指的平方根w吗？因为根据 的文档random，参数应该是“sigma”，即标准差。将 SD 定义为方差的平方根是很自然的，我猜它是w. 那么你在第一件作品中做得很好。

如果是这样，在你的第二段代码中，你为什么要sqrt()穿上w_sq？你的意思是取方差的第四根吗？从您对期望的定义来看，我认为这是不正确的。请看一下。

另一方面，是w_sq单个数字，还是数组？

  • 如果是数字——

tmpp / sqrt(w_sq(1))应该是tmpp / sqrt(w_sq)，尽管它们实际上并没有什么不同。

  • 如果它是一个数组 -

您可能希望将所有代码放在一个for循环中。循环遍历w_sq，每次它选择数组中的一个元素（将变量命名为 say w_sq_elem），然后让其余代码像单个数字一样执行操作。

无论如何，(log(x)-m)/sqrt(w_sq)并且tmpp / sqrt(w_sq(1))正在告诉不同的信息w_sq。第一个假设它是一个数字，所以除法可以简单地是 a /，而不是./。第二个表示它是一个数组，因此您正在选择它的第一个元素。但是数组在这里没有意义，因为根据我的理解，您不是将 10,000 个点除以 10,000 个x不同的方差。

m = 3;
w_sq = 2;
eta1 = -2;
ST = random('Lognormal',m,w_sq,10000000,1);
hlp = zeros(10000000,1);
hlp(ST<=2) = ST(ST<=2);
hlp(hlp>0) = hlp(hlp>0).^(eta1+1); % 0^(eta1+1) gives infinity
mean(hlp)

tmpp = integral(@(x) x.^(eta1) .* normpdf((log(x)-m)/w_sq),0,2) ;
tmpp / w_sq

>> untitled

ans =

    0.2944


ans =

    0.2948

>>
于 2014-07-17T02:00:14.187 回答