求助：weibull 数据生成

shiyiming

假定生存数据服从weibull分布（两参数a,b），已知中位生存时间为1年（或平均生存时间为2年），如何确定参数a,b，并生成生成数据。    t=rand('weibull',a,b);

shiyiming

approximately, a =0.733, b =1.648, based on mean =2 , median =1.
when sample size =10000, summary is as

[quote:2zuh5y9g]The MEANS Procedure

                                             Analysis Variable : x

                                                                                                                             
                          Mean             Median             Std Dev              Lower  Quartile               Upper Quartile
                  ----------------------------------------------------------------------------
                     1.9968633       0.9942942       2.7818681       0.2966301       2.5870673
                  ----------------------------------------------------------------------------[/quote:2zuh5y9g]

you can calculate a & b like

[code:2zuh5y9g]** supply initial values **;
data init;               
  a=.5; b=1;
run;

proc model data=init;
  eq.mean = b*gamma(1+1/a)-2;
  eq.median = b*(log(2))**(1/a) - 1;
  solve a b/out=solution;
run;

proc print data=solution;
run;[/code:2zuh5y9g]

JingJu

shiyiming

1.a,b的初始值是随意设定的吗？
2.下面是一篇文献的一段代码：
“假定生存时间服从威布尔分布，weibull(1,0.7);删失数据为右删失，删失分布为均匀分布U(0，c),删失比例设置为50%。产生样本量为100的生存时间的随机数据。通过估计，c=1.45。”
data lifetime;
a=1;b=0.7;c=1.45;
do i=1 to 1000;
censor=1;
t=rand('weibull',a,b);
h=rand('uniform')*c;
ht=min(t,h);
if (h < t) then censor=0;
output;
end;
run;
c是怎么确定的呢？c会随a,b,删失比例的变化而变化吧？
谢谢你！

shiyiming

应该是这样的。c好像要精确求解并不是很容易。但是显而易见，通过调整C来调节删失比例是很自然的事情。

shiyiming

jingju11大哥，我发现这样做出来的结果不是我想要的。没有删失数据的中位生存时间为1年，最后生成有删失的生存时间点后，中位生存时间就不是1年了，[color=#FF0000:1a6xmn2j]是不是应该根据终点时间来确定a,b.[/color:1a6xmn2j]
假定生存时间服从weibull分布，终点时间为1年，删失数据为均匀分布U（0，c），删失比例为50%，a,b怎样确定？

shiyiming

你无法准确知道生存的中位数，虽然你可以估计它，比如km。但是你利用某个分布去产生数据，其理想median是可知的。可不可以这么理解，如果你按照某个中位数xx的分布去产生数据，最后这个数据因为删失数据而导致估计的median偏离了原来的median。其实生存的中位数并未变化，只是因为估计方法的偏差而导致了不同的median。
京剧

shiyiming

[code:2n6zrajq]data lifetime;
a=0.733;b=1.648;c=3.5;
do i=1 to 1000;
censor=1;
t=rand('weibull',a,b);
h=rand('uniform')*c;
ht=min(t,h);
if (h < t) then censor=0;
output;
end;
run;[/code:2n6zrajq]
按上边的代码生成的数据出来的结果是：

[color=#FF0000:2n6zrajq]MEANS Procedure
变量   均值    中位数
t         2.045       1.067
ht       0.921        0.602[/color:2n6zrajq]

最后的生存时间应该是ht而不是t,因为我刚开始理解有误，生存时间不是正态分布，所以用均值估计参数是不对的，我想要的是最后ht的中位生存时间是1，结果显然不是，而且t值里有很大的值（大于15），是不是刚开始就应该设一个终点时间，比如说2年，所有的t值都应该是小于2的？

shiyiming

当然。这也是为什么要有survival 分析的原因：

proc lifetest;
time ht*censor(0);
run;

shiyiming

如果ht的中位生存时间为1年，a，b该怎么样确定啊？