任意の確率密度分布に従う乱数の生成（von Neumannの棄却法）

投稿日時: 2013年10月13日投稿者: pacocat

任意の分布に関して乱数を作りたかったのでメモ。
分布の逆関数を容易に求められる場合は、逆関数法を使ってしまうのが簡単なのだけど、

どんな（一価の）関数にも適用出来る
アルゴリズムが作りやすい

といった理由でここでは棄却法を使います。

以下のページを参考にしています。
数値計算法 2011/6/29 (pdf)
15 乱数とモンテカルロ法 (pdf)

von Neumannの棄却法

考え方はとてもシンプルで、

分布関数を含む範囲内で一様乱数のペア ( $x_{i}, y_{i}$ )を生成

具体的には分布関数の範囲を $[x_{min}, x_{max}], [0, y_{max}]$ 、
確率変数 $x_{rand}, y_{rand}$ を $0 \leq x_{rand}, y_{rand} \leq 1$ に従う一様分布と仮定した時に、

$x_{i} = x_{min} + (x_{max} - x_{min})x_{rand}\\ y_{i} = y_{max} y_{rand}\\$

を作る。

生成した乱数ペア( $x_{i}, y_{i}$ ) が分布関数内に入っていれば、 $x_{i}$ を分布に従う乱数として採用（下図でいうと赤点が採用、×が不採用）

アルゴリズム的にはこんな感じ？（コード中にあるテストファイルはこちらから→”random_dist.txt“）

結果のヒストグラム：

歪正規分布（Skew Normal Distribution）をmatplotlibでプロット

投稿日時: 2012年12月28日投稿者: pacocat

返信

正規分布やポワソン分布はいろんなパッケージが準備されてるけど、歪正規分布はなかなか見つからなかったので、
PASTEBIN
skew normal distribution in scipy
Skew Normal Distribution
などを参考にしてプロットしてみました。

歪正規分布（Skew Normal Distribution）

歪正規分布の確率密度関数はパラメータ $\lambda$ を使って以下のように定義されている。

$\displaystyle f(x,\lambda) = 2\phi (x) \Phi (\lambda x), \; \lambda >0$

ここで、

$\displaystyle \phi (x)=\frac{1}{\sqrt{2\pi}}\exp (-\frac{x^2}{2})$

$\displaystyle \Phi (x)=\int ^{x}_{-\infty} \phi (t)dt = \frac{1}{2}[1+erf (\frac{x}{\sqrt{2}})]$

コードはPASTEBINに紹介されているように定義されてる関数をそのまま返すだけ。erf関数はscipy.special.erfを使う。

例えばこんな感じ。

from pylab import plot,show
from numpy import exp,sqrt,pi,arange
from scipy.special import erf

# probability distribution function
def pdf(x):
    return exp(-x*x/2.0)/sqrt(2*pi)

# cumulative distribution function
def cdf(x):
    return (1.0+erf(x/sqrt(2)))/2.0

# skew normal distribution
def skew(x,mu,sigma,l):
    z= (x-mu)/sigma
    return 2/sigma * pdf(z) * cdf(l*x)

# define distribution parameters
mu=0.0
sigma=1.0
lam=4.0

# plot 
x=arange(-10,10,0.01)
y=skew(x,mu,sigma,lam)
plot(x,y)
show()

プロット

パラメータをいじってみたらこんな感じ。簡単のため平均=0、標準偏差=1にしてあります。

こっちは標準偏差を変えてみた場合（ $\lambda=3$ ）

超新星のSEDから光度曲線を計算（Example編）

投稿日時: 2012年11月27日投稿者: pacocat

返信

「超新星のSEDから光度曲線を計算」でNugent templateを使って計算した光度曲先例。
Bバンド最大光度は、以下を参考。
A Comparative Study of the Absolute Magnitude Distributions of Supernovae
Richardson et al. 2002（ADS）
# Habble定数が $H_0=60[{\rm km }{\rm s}^{-1} {\rm Mpc}^{-1}]$ になってるので注意。例えば $H_0=70$ に変換する時は、
$\displaystyle 5\log_{10}(\frac{70}{60})$ を加えること。

Ia型超新星（z=0.0, 0.5, 1.0）

Ibc型超新星（z=0.0）

IIL型超新星（z=0.0）

IIP型超新星（z=0.0）

IIn型超新星（z=0.0）

color evolution

超新星のSEDから光度曲線を計算

投稿日時: 2012年11月24日投稿者: pacocat

返信

計算式については「天体のfluxからBroad-band magnitudeを計算」を参照。

超新星のtemplate SED

Ia型はEric HsiaoのHPからダウンロード可能。
Peter NugentのHPからも、Type Ia 1991T-like, Type Ia 1991bg-like, Type Ib/c, Type Ib/c high-velocity, Type IIP, Type IIL, Type IInなどの各種templateが手に入る。

Filter Response

ここでは、Suprime-Cam Filtersをここからダウンロードして使っています。
# 使っているレスポンスは「フィルター + CCD + 主鏡反射率 + 主焦点補正光学系 + 大気透過率(secz=1.2) 」のもの。

計算手順

以下、HsiaoのType Ia template（snflux_1a.dat）をもとに光度曲線を計算する。

まずHsiaoのtemplateでは-20日から85日までのSEDが一つのファイルにまとまっているので、あらかじめspectrum_m20d.dat, spectrum_m19d.dat,…, spectrum_m01d.dat, spectrum_p00d.dat, spectrum_p01d.dat,…, spectrum_p85d.datのように各日毎にファイルを分離しておく。
後はABMAGの計算式、
$\displaystyle {\rm ABMAG}(R)=-2.5\log_{10}\frac{\int f_{\nu}\frac{R}{\nu}d\nu}{\int \frac{R}{\nu}d\nu}-48.6$
に従って計算するだけ。ただし、z>0の計算をする場合はスペクトルと時間が波長方向に(1+z)倍に伸びるので注意すること。

その他詳しい計算などは
K‐Corrections and Extinction Corrections for Type Ia Supernovae
Nugent et al. 2002 (ADS)
を参照。

コード例

distance modulusの計算にはPythonの宇宙論パッケージCosmoloPyを使用。
ここではVバンド絶対等級で-19.2のIa型超新星がSubaru/SCのB-,V-,Rc-,i’-,z’-bandでどのように観測されるかを計算している。

import numpy as np
import scipy.interpolate
import cosmolopy.distance as cosmocd

# convert wavelength to frequency
def convert_nu(filter):
    filter_lam=filter[:,0]*1.0e-8 # [cm]
    filter_nu=np.array([0.0] * len(filter_lam))
    for i in range(len(filter_lam)):
        filter_nu[i]=c/filter_lam[len(filter_lam)-1-i] # [Hz]
    return filter_nu*1.0e-15

# convert f_lambda to f_nu
def convert_fnu(filter):
    filter_lam=filter[:,0]*1.0e-8 # [cm]
    filter_responce=np.array([0.0] * len(filter_lam))
    for i in range(len(filter_lam)):
        filter_responce[i]=filter[len(filter_lam)-1-i,1]
    return filter_responce

# calculate ABmag
def abmagnitude(model_nu15,model_fnu,filter_nu15,filter_responce,z):
    flux=scipy.interpolate.interp1d(model_nu15,model_fnu,kind="linear")
    numerator=0.0
    denominator=0.0
    for i in range(len(filter_nu15)-1):
        for j in range(len(model_fnu)):
            if (filter_nu15[i]*(1.+z) = model_nu15[j-1]):
                dnu=(filter_nu15[i+1]-filter_nu15[i])
                numerator+=flux(filter_nu15[i]*(1.+z))*(1.+z)*filter_responce[i]/filter_nu15[i]*dnu
                denominator+=filter_responce[i]/filter_nu15[i]*dnu
    mag=-2.5*np.log10(numerator/denominator)-48.6
    return mag

# cosmological parameters
cosmo = {'omega_M_0':0.3, 'omega_lambda_0':0.7, 'omega_k_0':0.0, 'h':0.70}
# speed of light [cm s-1]
c = 2.989e10
# absolute V-band magnitude
MaV=-19.2

# load filter information
Bfilter=np.loadtxt('B-bandレスポンスファイル')
Vfilter=np.loadtxt('V-bandレスポンスファイル')
Rfilter=np.loadtxt('R-bandレスポンスファイル')
ifilter=np.loadtxt('i-bandレスポンスファイル')
zfilter=np.loadtxt('z-bandレスポンスファイル')

# define responce (as a function of nu)
Bfilter_nu15=convert_nu(Bfilter) # [E-15 Hz]
Bfilter_responce=convert_fnu(Bfilter)
Vfilter_nu15=convert_nu(Vfilter) # [E-15 Hz]
Vfilter_responce=convert_fnu(Vfilter)
Rfilter_nu15=convert_nu(Rfilter) # [E-15 Hz]
Rfilter_responce=convert_fnu(Rfilter)
ifilter_nu15=convert_nu(ifilter) # [E-15 Hz]
ifilter_responce=convert_fnu(ifilter)
zfilter_nu15=convert_nu(zfilter) # [E-15 Hz]
zfilter_responce=convert_fnu(zfilter)

# calculate AB magnitude from z=0 to z=1.2 in dz=0.05 interval
for z in np.arange(0.0,1.21,0.05):
    fout=open('z%03d' % (z*100.0),'w')
    if z>0:
        comoving_distance=cosmocd.comoving_distance(z, **cosmo) # comoving distance [Mpc]
        distance_modulus=5.*np.log10(comoving_distance)+5.*np.log10(1.+z)+25.
    else:
        comoving_distance=0.0
        distance_modulus=0.0
    # epoch loop
    for epoch in np.arange(-20,86,1):
        if epoch < 0:
            specdata='spectrum_m%02d.dat' % -epoch
        else:
            specdata='spectrum_p%02d.dat' % epoch
        model_lam=np.loadtxt(specdata)[:,0]*1.0e-8 # [cm]
        model_flam=np.loadtxt(specdata)[:,1]
        NBIN=len(model_lam)
        model_fnu= np.array([0.0] * NBIN)
        model_nu= np.array([0.0] * NBIN)
        for i in range(NBIN):
            model_fnu[i]=model_lam[NBIN-1-i]*model_lam[NBIN-1-i]*(model_flam[NBIN-1-i]*1.0e8)/c
            model_nu[i]=c/model_lam[NBIN-1-i]
        model_nu15=model_nu*1.0e-15
        # integration
        Bmag=abmagnitude(model_nu15,model_fnu,Bfilter_nu15,Bfilter_responce,z)+MaV+distance_modulus
        Vmag=abmagnitude(model_nu15,model_fnu,Vfilter_nu15,Vfilter_responce,z)+MaV+distance_modulus
        Rmag=abmagnitude(model_nu15,model_fnu,Rfilter_nu15,Rfilter_responce,z)+MaV+distance_modulus
        imag=abmagnitude(model_nu15,model_fnu,ifilter_nu15,ifilter_responce,z)+MaV+distance_modulus
        zmag=abmagnitude(model_nu15,model_fnu,zfilter_nu15,zfilter_responce,z)+MaV+distance_modulus
        print '%4.3f %5.1f %6.3f %6.3f %6.3f %6.3f %6.3f' % (z,epoch*(1.+z),Bmag,Vmag,Rmag,imag,zmag)

天体のfluxからBroad-band magnitudeを計算

投稿日時: 2012年11月24日投稿者: pacocat

返信

一度ちゃんとまとめておきたかったのでメモ。内容はLBLのNaoさんが院生時代に使っていたものを参考にしています。

まず、等級は天体のflux $f_{\lambda}$ , reference spectrum $f_{\lambda}^{*}$ , fileter response $R(\lambda)$ , reference magnitude $m_R^{*}$ を用いて以下のように定義できる。

$\displaystyle m_R=-2.5\log_{10}\frac{\int \lambda f_{\lambda}R(\lambda)d\lambda}{\int \lambda f_{\lambda}^{*}R(\lambda)d\lambda}+m_{R}^{*}$

Vega Mag, ST Mag, AB Magのうち、Vega MagはreferenceにVega fluxを使ったもの、ST MagとAB Magはreferenceとして一定の $f_{\lambda}, f_{\nu}$ を使ったものと定義されている。詳しくは以下。

VEGAmag : Magnitude system where Vega has magnitude 0 at all wavelengths by definition. The vega magnitude of a star with flux F is -2.5 log10 (F/F_vega) where F_vega is the current flux spectrum of Vega from the CALSPEC archive.
STmag and ABmag: Both systems define the absolute physical flux density for a point source. The conversion is chosen so that the magnitude at V corresponds roughly to that in the Johnson system. In the STmag system, the flux density is expressed per unit wavelength, while in the ABmag system, the flux density is expressed per unit frequency.
Hubble Space Telescope ACS Flux Calibration and Zeropointsから引用。

$f_{\lambda}$ based magnitude: ST Mag (Koorneef et al. 1986)

ST Magは $f_{\lambda}^{*}=3.63078 \times 10^{-9} [{\rm erg } {\rm cm }^{-2} s^{-1} {\rm \AA}^{-1}]$ で固定したもの。

このfluxでSTMAGが0になるようにするには、

$\displaystyle m_R^{*}=-2.5\log_{10}(f_{\lambda}^{*})=-21.10$

であればいいので、STMagは以下のように計算出来る。

$\displaystyle {\rm STMAG}(R)=-2.5\log_{10}\frac{\int \lambda f_{\lambda}Rd\lambda}{\int \lambda Rd\lambda}-21.10$

$f_{\nu}$ based magnitude: AB Mag (Oke 1964)

AB Magは $f_{\nu}^{*}=3.63078 \times 10^{-20} [{\rm erg } {\rm cm }^{-2} s^{-1} {\rm Hz}^{-1}]$ で固定したもの。

一番最初の定義式に戻って $\lambda$ から $\nu$ を以下のように変数変換する。

$\displaystyle \lambda=\frac{c}{\nu}$ $\displaystyle f_{\lambda}=f_{\nu}\frac{d\nu}{d\lambda}=f_{\nu}\frac{c}{\lambda^2}$

ABMAGが0になるようにするには、

$\displaystyle m_R^{*}=-2.5\log_{10}(f_{\nu}^{*})=-48.6$

であればいいので、結局ABMAGは以下のように計算出来る。

$\displaystyle {\rm ABMAG}(R)=-2.5\log_{10}\frac{\int f_{\nu}\frac{R}{\nu}d\nu}{\int \frac{R}{\nu}d\nu}-48.6$

また、ABMAG⇄STMAGのconversionは（上記の変数変換を使って）計算すると以下のようになる。

$\displaystyle {\rm ABMAG}(R)-{\rm STMAG}(R)=-2.5\log_{10}\frac{\int \lambda Rd\lambda}{\int \frac{R}{\lambda}d\lambda}+18.69$

これによって、ABMAGを $\lambda$ で表すことも可能。

$\displaystyle {\rm ABMAG}(R)={\rm STMAG}(R)-2.5\log_{10}\frac{\int \lambda Rd\lambda}{\int \frac{R}{\lambda}d\lambda}+18.69\\ =-2.5\log_{10}\frac{\int \lambda f_{\lambda}Rd\lambda}{\int \lambda Rd\lambda}-21.10-2.5\log_{10}\frac{\int \lambda Rd\lambda}{\int \frac{R}{\lambda}d\lambda}+18.69\\ =-2.5\log_{10}\frac{\int \lambda f_{\lambda}Rd\lambda}{\int \frac{R}{\lambda}d\lambda}-2.408$

Vega magnitude

最初の定義によりVega magnitudeは、

$\displaystyle {\rm Vega Mag}=-2.5\log_{10}\frac{\int \lambda f_{\lambda}Rd\lambda}{\int \lambda f_{\lambda}^{Vega}Rd\lambda}+0.0$

なので、

$\displaystyle {\rm ABMAG}-{\rm Vega Mag}=-2.5\log_{10}\frac{\int \lambda f_{\lambda}^{Vega}Rd\lambda}{\int \frac{R}{\lambda}d\lambda}-2.408$

従って、

$\displaystyle {\rm ABMAG}(R)={\rm Vega Mag}-2.5\log_{10}\frac{\int \lambda f_{\lambda}^{Vega}Rd\lambda}{\int \frac{R}{\lambda}d\lambda}-2.408$

wordpressにWP Latex、SyntaxHighlighter Evolvedを導入

投稿日時: 2012年11月24日投稿者: pacocat

返信

研究メモに使えそうなのでLatexプラグイン「WP Latex」と「SyntaxHighlighter Evolved」を入れてみた。

WP Latex

使い方は[ latex ]と[ /latex ]でくくった部分にtexコマンドを打ち込むだけ。以下のオプションが利用可能。

オプション	説明	例
color	色の変更	color=”000000″
background	背景色の変更
size	サイズの変更（可能なサイズは-4から4まで）0がnormal size（12pt）	size=”4″

例えば

[ latex size="4" ]\displaystyle \int ^{\infty}_{-\infty} e^{-ax^2}dx=\sqrt{\frac{\pi}{a}}[ /latex ]

と入力すれば以下のように表示される（入力の際[]内のスペースは無視）。
$\displaystyle \int ^{\infty}_{-\infty} e^{-ax^2}dx=\sqrt{\frac{\pi}{a}}$

その他詳しい説明はWP LaTeX FAQ参照。

SyntaxHighlighter Evolved

コードも貼りたいのでSyntaxHighlighter Evolvedもインストール。使い方は同様に[ python ][ /python ]などでソースコードを囲むだけ。

from itertools import count
from math import sqrt

def prime_gen():
    primes = []
    for n in count(2):
        if all(n%p for p in primes if p <= sqrt(n)):
            primes.append(n)
            yield n

g = prime_gen()
for i in xrange(1000):
    print g.next()

References:
ソースコードをキレイに表示するWordPressプラグイン「SyntaxHighlighter Evolved」
Clojure prime numbers lazy sequence

Pacocat's Life

Every day is a new day

カテゴリーアーカイブ: 研究

任意の確率密度分布に従う乱数の生成（von Neumannの棄却法）

von Neumannの棄却法

歪正規分布（Skew Normal Distribution）をmatplotlibでプロット

歪正規分布（Skew Normal Distribution）

プロット

超新星のSEDから光度曲線を計算（Example編）

Ia型超新星（z=0.0, 0.5, 1.0）

Ibc型超新星（z=0.0）

IIL型超新星（z=0.0）

IIP型超新星（z=0.0）

IIn型超新星（z=0.0）

color evolution

超新星のSEDから光度曲線を計算

超新星のtemplate SED

Filter Response

計算手順

コード例

天体のfluxからBroad-band magnitudeを計算

$f_{\lambda}$ based magnitude: ST Mag (Koorneef et al. 1986)

$f_{\nu}$ based magnitude: AB Mag (Oke 1964)

Vega magnitude

wordpressにWP Latex、SyntaxHighlighter Evolvedを導入

WP Latex

SyntaxHighlighter Evolved