给定任意一个概率密度函数,生成随机数

工作中经常遇到这样一个问题:给定一个概率密度分布函数pdf,要生成N的随机数满足这个pdf。比如说生成一个线性的pdf,较大的x出现的概率大一些,较小的x出现的机会小一些,线性比例分配。MATLAB不提供这样的函数。

我写了一个简单的程序,利用Gibbs sampler的办法撒点得到,效率还不错。对于一个线性分布来说,生成10000个点需要1.9秒,10万个点仅需要不到19秒的时间。测试的机器是Dell 4核没做并行计算。

代码贴在这里,也提供文件下载frnd.m


% You really need FRND(x,p,a,b,N) if you want to generate random numbers
% following someone known possibility distribution function.
%
% p is the PDF of x at the range of [a,b]
% FRND(x,p,a,b,N) retures N random numbers following distribution of PDF
% x and p should be column vectors with the same length.
%
% x: a column vector
% p: a vector with the same size of x. The p is the possibility at value x.
% a: lower bound of random numbers, a scalar
% b: upper bound of random numbers, a scalar
% N: total count of result you want, a scalar
% result: a column vector with length N, which follows PDF defined by (x,p)
%
% by Shuang Gao
% please contact sgao@nao.cas.cn if you have any questions.
function xf = frnd(x,p,a,b,N)
maxp = max(p);
minp = min(p);

xf = zeros(N,1);
for i = 1:N;

xtry = unifrnd(a,b);
ptry = unifrnd(minp,maxp);
y = interp1q(x,p,xtry);
while ptry>y,
xtry = unifrnd(a,b);
ptry = unifrnd(minp,maxp);
y = interp1q(x,p,xtry);
end
xf(i) = xtry;
end

用太阳临近星团限制恒星形成率

最近的Oberseminar要讲一下这个,Constraining the star formation rate in the Solar neighborhood with star cluster arXiv:1104.2182v1 (文章中涉及到星团和质量流失的描述很简单,本人通过文献补充了一些基本背景。)
下载ppt (pdf文件,其实不能叫ppt)

背景介绍

恒星星团是星系中靠引力束缚着的恒星集团,依据目前的观测,可以分为以下几类:球状星团GC,疏散星团OC,年轻大质量星团YMC,嵌埋星团(embedded cluster EC)。

  • GC包含百万量级的低质量老年恒星,存在于晕中,自引力强,呈球状;
  • OC包含几百到几万量级的不同类型恒星,存在于盘中,形状不规则,自引力弱;
  • YMC包含几百到几万量级的年轻大质量恒星,存在于盘中,恒星平均质量>5 太阳质量;
  • EC是OC的早期阶段,诞生于巨分子云(GMC)中,光学难以探测。

恒星形成的历史,用单位时间形成的恒星质量来描述,它是时间的函数SFR(t)。太阳临近空间的本地SFR对于构建银河系模型、星族研究、了解银河系的历史等都非常有用处。本文是用星团的年龄分布限制SFR。

星团诞生后,会经过三个演化阶段:嵌埋阶段、恒星演化阶段、动力学演化阶段。Lada & Lada (2006)提出所有的星团都产生于GMC,也就是说星团诞生时都是嵌埋星团,但大部分(90%以上)在这个阶段就死亡了。第二个阶段是星团中的成员星独自演化的阶段,演化过程中星风带走质量,演化结果是白矮星、中子星或黑洞。第三个阶段和第二阶段并没有明确界限,星团受到几种动力学干扰也会走向死亡。

造成星团死亡(解散 dissolution)的原因就是质量丢失,质量丢失一共有6种机制:

  • 恒星演化:通过星风带走质量,白矮星、中子星、黑洞部分被踢出;
  • 银河系的潮汐力瓦解:变星,拉扯出部分恒星(形成星流);
  • 旋臂产生的激波:激波作用时间越长效果越严重;
  • 遭遇巨分子云(GMC):包括嵌埋星团等;
  • 抛射(ejection):单个恒星近距离交互时的N body作用获得的额外速度大于逃逸速度后逃出;
  • 蒸发(evaporation):大量恒星充分交互后发生质量分层(mass segmetation),恒星获得相近的速度,大质量星速度弥散小,小质量星速度弥散大向外层迁移最终逃出。

针对这6种质量丢失的机制,有大量文献给出了观测的和N体数值模拟的估计,它们各自有不同的质量丢失率dM/dt。总的星团质量丢失率是这6种质量丢失率简单相加。如果我们把条件限定在太阳附近,就可以给出完整的6种质量丢失率的表达式对时间的函数。用一个Runge-Kutta方法解方程就可以得到M随着t的变化M(t)。

观测与模拟ADF

在距离太阳1kpc的范围内共计442个星团,它们的年龄分布可以完整地得到。这是观测ADF。与之对比的是一个模拟的ADF。模拟方法如下:

  1. 按照某个给定的初始质量函数随机分配质量,随机分配年龄ta;
  2. 根据上面计算的M(t)计算每个星团演化到自己年龄ta时的值M(ta);
  3. 重复这样的步骤直到总质量满足预设的一个SFR对时间的积分;
  4. 演化后的M(ta)小于100太阳质量的忽略不计,对大于100太阳质量的星团计数,得到模拟ADF。
  5. 上面的步骤重复100次平均用以降低随机误差。

如上图,横轴是百万年的年龄的对数尺度,数轴是单位时间间隔的星团数。带误差棒的原点是观测ADF,其他线都是模拟ADF。其中

  • 最粗的点划线是假设了常数的SFR,在两个时间范围内与观测偏差很大,一个是200到600万年,1个是最近;
  • 实线是符合最好的SFR,在这两个时间范围内调高了SFR的值,具体数值对应在图下面的子图中;
  • 两外两条线是这个SFR上下20%偏差的位置;阴影是模拟100次的结果的偏差范围。

结论与讨论

这个最好的SFR的结果告诉了我们以下几件事

  • 在最近(9百万年之内)和200到600万年之前,本地恒星形成有较高的速率;
  • 平均SFR是2500太阳质量每百万年
  • 常数的SFR不能解释本地星团

在模拟中发现,共用了330000个本地星团,但是

  • 91.2%都在最初10Myr就毁灭了;
  • 这个数字与Lada等人2006年用嵌埋星团计算的结果一致。
  • 这个结论暗示了,最初10Myr内毁灭的全部是嵌埋星团。

延伸阅读

http://en.wikipedia.org/wiki/Star_cluster

http://en.wikipedia.org/wiki/Open_cluster

Embedded Cluster review

Young massive cluster review

Levenberg-Marquardt方法

Levenberg-Marquardt方法是Levenberg和Marquardt两个人分别在1944年和1963年提出的优化算法。这种方法专门用于计算函数F(x)的最小值。适用的F(x)是一系列函数f(x)的平方和的形式。

F(x)=frac{1}{2}sum_{i=1}^m{[f(x)]^2}

显然,这是一个非线性函数。

这种寻找最线性函数“最小值”的方法,很自然地可以应用在参数拟合上。假如我们定义函数与数据之间的chi^2值如下:

chi^2=frac{(ln{data}-ln{model})^2}{sigma^2}

这显然是一个平方项。如果函数是一系列chi^2值的和,就恰好可以应用Levenberg-Marquardt方法做参数拟合。例如若干个星族之和构成的星族合成模型与观测的比较,每个组分分别考虑为一个单独的chi^2值,总的结果可以通过Levenberg-Marquardt方法拟合。针对天文的观测数据,data往往作为拟合的参照标准,拟合的权重sigma^2通常定义为data的倒数值。

Fortran程序在NR中给出了完整的子程序mrqmin和mrqcof。前者计算拟合的参数值即最小误差值所在的位置,后者计算误差值本身。

具体程序在这里:http://hahn.lbl.gov/svn/rsxap/trunk/source/NR/mrqmin.f

在matlab的优化工具箱里,提供了5种最小卡方值的拟合算法。Levenberg-Marquardt方法是其中之一。在matlab中运行optimtool即可打开优化工具箱,其中可以选择非线性拟合的方法为Levenberg-Marquardt。

论文中的语言问题

按照老板给的修改,总结一下常出错的几个地方。

the不能缺少的地方:

the Einasto laws (以某人命名的定律),the SDSS ,the g band

不能加the的地方:

TRILEGAL model,de-reddened magnitude

其他习惯用语:

instead 作为副词 be used instead.
instead of :
They burned instead of burying their dead. 而不是
Instead of money he gave promises. 相当于for
Q is a Letter of small use, and put only in stead of C. 替换某某位置

combination代替pair形容SFR and AVR

这两个小东西后面要有逗号:
..., e.g., ...
..., i.e., ...

不用章节做主语:
In Section 6 we summarise ...而不是Section 6 summarises

介词:
for each direction
on the sky
in the direction of the sky,in the NGP direction

upper right 和lower left (不是top right 和 bottom left)

as described in Section 2

refered to. (to不能省略)

 

其他:

specified不是special

best fit 或者 best-fitting

所有的mean都换成stand for