请各位高手帮我看看程序

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

1. %let m=200;
   
2. %let n=200;
   
3. %let gamm=2;
   
4. %let R0=0.9;
   
5. %let lambda=0.5;
   
6. %let seed=1;
   
7. %let Fdis=FINV(1-&R0,2,2*&gamm);
   
8. %let rho=&gamm/&Fdis;
   
9. %let tau=1.1;
   
10. 
    
11. 
    
12. proc iml;
    
13. 
    
14. reset noprint;
    
15. 
    
16. 
    
17. randx=ranexp(repeat(&seed,&m))/λ
    
18. 
    
19. randy=&tau*rangam(repeat(&seed,&n),&gamm);
    
20. 
    
21. one1=j(&m,1,1);
    
22. u=randx`*one1;
    
23. one2=j(&n,1,1);
    
24. v=randy`*one2;
    
25. 
    
26. 
    
27. start inner(t) global(x);
    
28. _vu=v/u; _uv =u/v; _gn =(&n-1)*&gamm;
    
29. _t =t*_vu;
    
30. B= beta(&gamm,(&n-1)*&gamm);
    
31. z =1/B *((_uv*_t) **(&gamm -1) *(1-(_uv*_t)) **(_gn -1)) *_uv *(1-x) **(&m -2) *(&m -1);
    
32. return(z);
    
33. finish;
    
34. start outer(s) global(x);
    
35. x =s;
    
36. t_interval = 0 ||x;
    
37. call quad(w, "inner", t_interval);
    
38. return (w);
    
39. finish;
    
40. /* when u <=v*/
    
41. if u < v then
    
42. do;
    
43. s1_interval= {0 1};
    
44. call quad(R0, "outer", s1_interval);
    
45. print R0;
    
46. end;
    
47. 
    
48. if u > v then
    
49. do;
    
50. _vu=v/u;
    
51. s2_interval= 0 ||_vu;
    
52. call quad(R1, "outer", s2_interval);
    
53. R2 =R1 -(1-_vu)**(&m-1);
    
54. print R2;
    
55. end;
    
56. 
    
57. 
    
58. quit;
    
59. %put _user_;

复制代码

运行不出结果，望各位赐教。谢谢！

扫码加我 拉你入群

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

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

分享0 收藏0 回帖

关键词：Interval RETURN Global finish repeat repeat inner start 程序

完全无懂…………

1# saiyasaibing

1. 
   
2. %let m=200;
   
3. %let n=200;
   
4. %let gamm=2;
   
5. %let R0=0.9;
   
6. %let lambda=0.5;
   
7. %let seed=1;
   
8. %let Fdis=%sysfunc(FINV(1-&R0, 2, 2*&gamm));
   
9. %let rho=%sysevalf(&gamm /&Fdis); %put ρ
   
10. %let tau=1.1;
    
11. proc iml;
    
12. reset noprint;
    
13. randx =ranexp(repeat(&seed,&m))/λ
    
14. randy =&tau *rangam(repeat(&seed,&n),&gamm);
    
15. one1=j(&m,1,1);
    
16. u =randx`*one1;
    
17. one2 =j(&n,1,1);
    
18. v =randy`*one2;
    
19. B =beta(&gamm, (&n-1)*&gamm);
    
20. _vu =v/u;
    
21. _gn =(&n-1)*&gamm;
    
22. print u v b _gn;
    
23. 
    
24. start inner(t) global(x, v, u, b, _vu, _gn);  
    
25.   z =1/b*(t **(&gamm -1) *(1-t) **(_gn -1)) *_vu *(1-_vu*x) **(&m -2) *(&m -1);
    
26.   return(z);
    
27. finish;
    
28. start outer(s) global(x);
    
29.   x =s;
    
30.   t_interval = 0 ||x;
    
31.   call quad(w, "inner", t_interval) eps=1E-10 peak=0.0001;  
    
32.   return (w);
    
33. finish;
    
34. if u <=v then do;
    
35.   s1_interval= 0 ||_vu;
    
36.   call quad(R0, "outer", s1_interval) eps=1E-10 peak=0.0001;
    
37. end; print u v R0[format=E21.15];
    
38. /* u =450; _vu =v/u;*/
    
39. if u > v then do;
    
40.   s2_interval= {0 1};
    
41.   call quad(R1, "outer", s2_interval) eps=1E-10 peak=1e-7;
    
42.   R2 =R1 -(1-_vu)**(&m-1);
    
43. end; print u v R2[format=E21.15];
    
44. quit;

复制代码

I modify the formula since I was using a false transformation before (I forgot the basic calculus even).
Even though I make the code execute, but, honestly, I am not clear what happens for SAS calculation the integral behind the scene, nor the usage of eps, peak and so on. I would ask someone who really knows about that. On the other hand, probably you can ask for a check on MatLab since so many people have mentioned its superiority on integration.
JingJu

3# jingju11
Thank you jingju, I definitely will try Matlab. I just don't wanna give up on SAS.