SAS可以直接产生指定偏度，峰度的分布数据吗？

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

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

立即领取

感谢您参与论坛问题回答

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

+2 论坛币

如题： SAS可以直接产生指定偏度，峰度的分布数据吗？

如：想办法产生，mean= 2, std=3, kurtosis=4,skewness=5 的这样一个分布的数据吗？

扫码加我 拉你入群

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

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

分享0 收藏0 回帖

关键词：Kurtosis skewness mean Skew RTO

还没有软件有这个功能吧？ 信息太少了

如果偏度，峰度均为0时，是正态分布。那么此时，是可以用rand，rannor等函数实现的。
如查不为0，我就不清楚是否有办法模拟了。

标准正态分布只有两个参数便可固定，即均值和方差
但是非标准的，用均值、方差，再加上偏度和峰度则不然，同样的偏度和峰度，可以有多个分布
如果我没有记错的话，应该是这样
具体还需要什么参数，我就不清楚了

刚我搜了下，下面两个链接可以看一下

https://stat.ethz.ch/pipermail/r-sig-finance/2008q3/002825.html

http://stackoverflow.com/questio ... w-and-kurtosis-in-r

Referring to Wicklin's "simulating data with sas"  and the Fleishman's paper  about cubic transformation method in Psychometrika(1978), i am mainly programing  in FCMP, with the help pf Macro facility.
The partial code in FCMP is as below.
In this case, my macro can be called as: Simu_Skew_Kurt(,2 3 4 5,,,100);
let me know if anything wrong.
More Details can be seen in my blog at  http://blog.sina.com.cn/s/blog_a3a926360101f6us.html

1. proc fcmp;
   
2. 
   
3.       %* input coefficients and output targeted moments: variance, skewness, kurtosis;
   
4.       subroutine Fleishman(c[*], x[*]) varargs; outargs x;
   
5.          x[1] = (c[1]**2) + 6*(c[1]*c[3]) + 2*(c[2]**2) + 15*(c[3]**2); /* variance */
   
6.          x[2] = 2*c[2]*((c[1]**2) + 24*(c[1]*c[3]) + 105*(c[3]**2) + 2); /* skewness */
   
7.          x[3] = 24*((c[1]*c[3])+(c[2]**2)*(1+(c[1]**2)+28*(c[1]*c[3]))+(c[3]**2)*(12+48*(c[1]*c[3])+141*(c[2]**2)+225*(c[3]**2))); /* excess kurtosis */
   
8.       endsub;
   
9. 
   
10.    %* input coefficient and output J matrix of Jacobian;
    
11.    subroutine FlDeriv(x[*],J[3,3]) varargs; outargs J;
    
12.       J[1,1] = 2*x[1] + 6*x[3];
    
13.       J[1,2] = 4*x[2];
    
14.       J[1,3] = 6*x[1] + 30*x[3];
    
15.       J[2,1] = 4*x[2]*(x[1] + 12*x[3]);
    
16.       J[2,2] = 2*((x[1]**2) + 24*(x[1]*x[3]) + 105*(x[3]**2) + 2);
    
17.       J[2,3] = 4*x[2]*(12*x[1] + 105*x[3]);
    
18.       J[3,1] = 24*(x[3] + (x[2]**2)*(2*x[1] + 28*x[3]) + 48*x[3]**3);
    
19.       J[3,2] = 48*x[2]*(1 + (x[1]**2) + 28*(x[1]*x[3]) + 141*(x[3]**2));
    
20.       J[3,3] = 24*(x[1]+28*x[1]*(x[2]**2)+2*x[3]*(12+48*(x[1]*x[3])+141*(x[2]**2)+225*(x[3]**2))+(x[3]**2)*(48*x[1]+450*x[3]));
    
21.    endsub;
    
22. 
    
23.    %* Iterations of Newton algorithm to resolve cubic coefficients;
    
24.    array cf[1,4]/nosymbols;
    
25.    rc =read_array('work.moments4',cf);
    
26.    array tr[3];
    
27.    %* Assign Targeted moments parameters;
    
28.    array ti[3]  (1 . .); ti[2] =cf[1,3]; ti[3] =cf[1,4];
    
29.    %* Fleishman coefficients with an initial guess by a polynomial regression;
    
30.    array x[3];
    
31.    x[1] = 0.95357 -0.05679*ti[3] +0.03520*ti[2]**2 +0.00133*ti[3]**2 ;
    
32.    x[2] = 0.10007*ti[2] +0.00844*ti[2]**3                            ;
    
33.    x[3] = 0.30978 -0.31655*x[1]                                      ;
    
34.    array J[3,3]   ;  %*Jacobian matrix;
    
35.    array re[3,3]  ;  %*Intermediate results;   
    
36.    array f[3]     ;  %*Objective function;
    
37.    array ep[3]    ;  %*Error updates;
    
38. 
    
39.    maxIter =&maxIter.; absConv =&absConv.;
    
40.    %* Initialize f- objective function;
    
41.    call Fleishman(x,tr);
    
42.    call subtractMatrix(tr, ti, f);
    
43. 
    
44.    mf =absConv +1;
    
45.    do iter =1 to maxIter by 1 while(mf >absConv);
    
46.       %* Solve error function;
    
47.       call FlDeriv(x, J);      
    
48.       call inv(J, re);
    
49.       call mult(re,f, ep);
    
50.       call mult(ep, -1, ep);
    
51.       %* Updating coefficients x;
    
52.       call addMatrix(x, ep, x);
    
53.       %* Updating objection fucntion f;
    
54.       call Fleishman(x,tr);
    
55.       call SubtractMatrix(tr,ti, f);
    
56.       %* Check error in Error allowance;
    
57.       mf =0;
    
58.       do i =1 to 3;         
    
59.          if abs(f[i])>mf then mf =abs(f[i]);
    
60.          end;
    
61.      end;
    
62.    %* Set status variable for if converging;
    
63.    call symputx('status', catx(' ', 'Converge in iterations of', iter-1));
    
64.    if mf > absConv then do;
    
65.       put "!!!!NOT converge in maximum iterations of" maxIter;
    
66.       call symputx('status', catx(' ', 'Not Coverge in iterations of', iter-1));
    
67.       end;
    
68.    %* output solved Fleiman coefficients;
    
69.    rc =write_array('work.Fleishman_cubic_coeffs',x);
    
70.    run;

复制代码

Jingju

jingju11 发表于 2013-8-23 05:31
Referring to Wicklin's "simulating data with sas"  and the Fleishman's paper  about cubic transforma ...

京剧大侠，问下你的方法是用方差、偏度和峰度这三个不同阶数的矩列方程
然后用Jacobian求解吗？
那么解就不是唯一的了，就是说分布都不是唯一的
然后在这个随机的分布上再随机得到一组基于该分布的观测值吗？
谢谢！

jingju11 发表于 2013-8-23 05:31
Referring to Wicklin's "simulating data with sas"  and the Fleishman's paper  about cubic transforma ...

非常感谢京剧大哥。我得好好消化一下。

I am not fully understanding what Fleishman's cubic transformation? I would think it is an approximation. On the other hand, we need 3 coefficients and we have 3 functions. Looks like the solution is unique somehow. JingJu

个人觉得分布函数是决定唯一分布的，楼上也提到了，指定了均值、方差、峰度、偏度的分布并不唯一。

我的想法是从一个分布出发，适当修改分布函数，实现目的，不过工作只完成了一半，偏度的修改还没有一个很好的方法。

从正态分布出发，其峰度的计算如下：

为了简便，这里只计算N(0,1)，N(mu,sig)的计算相似。

Ps:按照定义计算，正态分布的峰值为3，很多软件都作了减3的处理。

修改密度函数，并计算均值、方差，

若给定峰度Kurt、均值Mu、方差Sig2，生成随机数

计算出:

k=Kurt/3;

mu=Mu/k;

sig2=Sig2/k;

转化为：

以密度函数

参数分别为(k,mu,sig2)生成随机数。

偏度的推导：

还没想到合适的修改方法。

这里只是给出一个想法，有兴趣的朋友可以研究一下，不局限于正态分布，我觉得这是个很有意思的问题。