请教关于r 和统计的 问题 问很多人都不会

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

1. > # the likelihood for parameter sigma and rho
   
2. > bvn_likelihood = function(x,sigma,rho){
   
3. +  n = length(x)/2; a = x[,1]; b = x[,2]
   
4. +        p = matrix(0,1,n)
   
5. +  for (i in 1:n){
   
6. +        ab =(a^2 - 2*rho * a * b + b^2)/(2 * (1 - rho^2)*sigma^2)
   
7. +        p = exp(-ab)/(2*pi*sigma^2*sqrt(1-sigma^2))
   
8. +       }
   
9. + prod(p)
   
10. + }
    
11. >  
    
12. >  
    
13. > #the prior density for sigma and rho
    
14. > bvn_prior = function(sigma, rho)  exp(-sigma)/(4*sqrt(abs(rho)))
    
15. >
    
16. >
    
17. > #the normalizing constant for the posterior distribution
    
18. > bvn_normalize = function(x){
    
19. + integrate2 = function(f, sigma.lower, sigma.upper, rho.lower, rho.upper) {
    
20. +  f2 = function(rho) integrate(f, lower = sigma.lower, upper = sigma.upper, rho = rho)$value
    
21. +  f3 = function(rho) sapply(rho, f2)
    
22. +  integrate(f3, rho.lower, rho.upper)
    
23. +  }
    
24. +  f = function(sigma, rho) bvn_likelihood(x=x,sigma, rho)*bvn_prior(sigma, rho)
    
25. +  integrate2(f, 0, Inf, 0, 1)$value +integrate2(f, 0, Inf, -1, 0)$value
    
26. +  a = integrate2(f, 0, Inf, 0, 1)$value; b = integrate2(f, 0, Inf, -1, 0)$value
    
27. +  a+b
    
28. + }
    
29. >
    
30. > #the marginal posterior density of rho
    
31. > bvn_posterior_rho = function(x, rho){
    
32. +  integrate (function(sigma) {
    
33. +    bvn_normalize(x)*bvn_likelihood(x, sigma, rho)*bvn_prior(sigma, rho)}, 0, Inf) $value
    
34. + }
    
35. >
    
36. > # the maginal prior density for rho
    
37. > bvn_prior_rho = function(rho)  integrate(function(sigma) bvn_prior(sigma, rho), 0, Inf)
    
38. >
    
39. > set.seed(1)
    
40. > n <- 40
    
41. > x1 <- rnorm(n); x2 <- rnorm(n); z <- rnorm(n)
    
42. > Xa <- cbind (x1, x2); Xb <- cbind (x1+0.5*z, x2+0.5*z)
    
43. > sigma_a = sd(Xa); sigma_b = sd(Xb)
    
44. >
    
45. > bvn_normalize(Xa)

复制代码

Error in integrate(f, lower = sigma.lower, upper = sigma.upper, rho = rho) :
  non-finite function value
In addition: There were 50 or more warnings (use warnings() to see the first 50)
请问 有哪位大神会啊～～知道哪里错了 能帮忙改一下 问了好多人了 都不会啊！！！

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：distribution Likelihood integrate Parameter function 统计

蓝色bao风雪 发表于 2014-11-12 22:52
&gt; # the likelihood for parameter sigma and rho
&gt; bvn_likelihood = function(x,sigma,rho){
+  n = le ...

integrate (function(sigma) {
bvn_normalize(x)*bvn_likelihood(x, sigma, rho)*bvn_prior(sigma, rho)}, 0, Inf)
这里感觉不对，其他的我不知道只看到了这里，我觉得你可以把你的思路，目标简单的描述下，而不只是代码。

ribbon 发表于 2014-11-12 23:34
integrate (function(sigma) {
bvn_normalize(x)*bvn_likelihood(x, sigma, rho)*bvn_prior(sigma, rho) ...

目的是要求bayes 的常数 第bvn_likelihood 是关于sigama和rho的似然函数 X 是二元正态分布 mean ＝ （0，0）
sigma 是x[1]和x［2］的方差 （都是sigma） 未知， correlation rho 也是未知
x~(X11;X12); (X21;X22); : : : ; (Xn1;Xn2). ~Normal所以，

第二个bvn_prior是求关于sigma和rho的先验密度
第三个bvn_normalize 是求常数c 等于 关于sigma 和rho的 似然函数＊先验密度 二次积分
就是这个地方好像出问题了

ribbon 发表于 2014-11-12 23:34
integrate (function(sigma) {
bvn_normalize(x)*bvn_likelihood(x, sigma, rho)*bvn_prior(sigma, rho) ...

> #the marginal posterior density of rho
> bvn_posterior_rho = function(x, rho){
+  integrate (function(sigma) {
+    bvn_normalize(x)*bvn_likelihood(x, sigma, rho)*bvn_prior(sigma, rho)}, 0, Inf) $value
+ }
>> bvn_posterior_rho(Xa,0.5)
Error in integrate(f, lower = sigma.lower, upper = sigma.upper, rho = rho) :
  non-finite function value
In addition: There were 50 or more warnings (use warnings() to see the first 50)
这个也不对  这个是求关于rho的后验密度

要不晚上一起讨论,如果你有时间.
把你的联系方式发给我.

请问最后解决了么？