請教關(guān)于r 和統(tǒng)計的 問題 問很多人都不會

是 否 +2 論壇幣

k人 參與回答

經(jīng)管之家送您一份

應(yīng)屆畢業(yè)生專屬福利!

求職就業(yè)群

趙安豆老師微信：zhaoandou666

經(jīng)管之家聯(lián)合CDA

送您一個全額獎學(xué)金名額~ !

立即領(lǐng)取

感謝您參與論壇問題回答

經(jīng)管之家送您兩個論壇幣！

+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)

復(fù)制代碼

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 回帖

關(guān)鍵詞：distribution Likelihood integrate Parameter function 統(tǒng)計

藍色bao風雪 發(fā)表于 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 發(fā)表于 2014-11-12 23:34
integrate (function(sigma) {
bvn_normalize(x)*bvn_likelihood(x, sigma, rho)*bvn_prior(sigma, rho) ...

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

第二個bvn_prior是求關(guān)于sigma和rho的先驗密度
第三個bvn_normalize 是求常數(shù)c 等于 關(guān)于sigma 和rho的 似然函數(shù)＊先驗密度 二次積分
就是這個地方好像出問題了

ribbon 發(fā)表于 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)
這個也不對  這個是求關(guān)于rho的后驗密度

要不晚上一起討論,如果你有時間.
把你的聯(lián)系方式發(fā)給我.

請問最后解決了么？