精华内容
下载资源
问答
  • 在上一篇文章(参见:使用模板库进行参数估计)中,我们使用网格搜索来找两个啁啾信号参数中最有可能的值,效果很好。但是如果我们想一起估算函数的所有参数怎么办?也许用更密集和/或更宽的网格值可以做到这一点,...

    70dbec6b863ccf1c87a09e0cf2e483fa.png

    在上一篇文章(参见:使用模板库进行参数估计)中,我们使用网格搜索来找两个啁啾信号参数中最有可能的值,效果很好。但是如果我们想一起估算函数的所有参数怎么办?也许用更密集和/或更宽的网格值可以做到这一点,但是计算很快就会成为问题。

    我们需要一种更聪明的方法来对参数的似然函数进行抽样--蒙特卡洛技术。我们将使用马尔可夫链蒙特卡洛(Markov Chain Monte Carlo, MCMC)方法(参见文末附录)执行参数估计,而pymc软件包可方便地提供此功能。

    MCMC就是从概率分布中选择随机样本,但选择方式应使最终的样本集近似于该概率分布。这种方法的优点是它能快速地找到概率分布的最大值,并且样本大多集中在这个最大值附近。其余的参数空间被忽略了,但我们通常并不关心。下面将展示如何在Python中执行MCMC采样以进行参数估计。

    01

    生成信号与噪声

    首先产生啁啾信号,然后模拟噪声,并将信号注入其中。

     1import numpy as np
    2import pymc
    3import pandas as pd
    4import seaborn as sns
    5from matplotlib import pyplot as plt
    6%matplotlib inline
    7np.random.seed(0# For reproducibility
    8
    9# We will use the same chirp function as previously
    10def chirp(t, tc, offset, A, dA, f, df):
    11    """12    Generate a chirp signal.    13    Arguments:14    t      -- array-like containing the times at which to evaluate the signal.15    tc     -- time of coalescence, after which the signal is terminated; 16              for times beyond tc the signal is set to the constant background.17    offset -- constant background term18    A      -- initial signal amplitude at t=019    dA     -- linear coefficient describing the increase of the amplitude with time20    f      -- initial signal frequency at t=021    df     -- linear coefficient describing the increase of the frequency with time22    """
    23    chirp = offset + (A + dA*t) * np.sin(2*np.pi*(f + df*t) * t)
    24    chirp[t>tc] = offset
    25    return chirp
    26
    27# Let's choose some values for our injected signal
    28tc_true     = 75
    29offset_true = 30
    30A_true      = 6
    31dA_true     = 0.05
    32f_true      = 0.2
    33df_true     = 0.007
    34
    35# Noise strength; keep it low here so we don't have to run sampling for too long
    36sigma = 50
    37
    38# Time axis
    39t = np.linspace(0,100,10001)
    40
    41# Injecting our signal into the noise
    42y_true = chirp(t, tc_true, offset_true, A_true, dA_true, f_true, df_true)
    43y_obs = np.random.normal(y_true, sigma)

    可视化信号与观测数据。

    1# Let's have a look
    2plt.plot(t, y_obs, color='y')
    3plt.plot(t, y_true, color='r');
    4plt.xlabel("t")
    5plt.ylabel("A")
    6plt.show()

    b3f3986b1a3dec2d951529c82752aeaa.png

    为了简化,与上篇文章相比,这次的噪音没有那么大,也没有使用较小的信号,因为我们不想为MCMC算法运行而等待太久。

    02

    设定pymc

    PyMC是一个可以实现多种贝叶斯统计方法的程序包,它提供了一些可以拟合概率模型的方法,最主要的拟合模型方法是MCMC算法。它要求你通过定义所有要估计的参数来构建模型,随后它会处理所有繁重的工作。

    我们需要将所有参数指定为独立的随机变量。对于每个参数指定:

    1. 先验分配的类型,即我们开始之前所知道的信息。在这种情况下,我们假设没有先验信息,因此我们只使用均匀分布。在实际应用中,使用先验信息,例如黑洞或中子星在宇宙中的分布可能会带来好处。

    2. 参数的名称,用于以后引用它。

    3. 均匀分布的取值范围。

    4. 初始值,我们选择为随机值。

    首先定义各参数为均匀分布

     1# Helper data to keep track of everything
    2parameters = ['tc''offset''A''dA''f''df']
    3bounds = {'tc':(0,100), 'offset':(0,100), 'A':(0,10), 'dA':(0,0.1), 'f':(0,1), 'df':(0,0.1)}
    4
    5# Defining our parameters using pymc.Uniform()
    6offset = pymc.Uniform('offset', bounds['offset'][0], bounds['offset'][1], \
    7         value = bounds['offset'][0] + (bounds['offset'][1]-bounds['offset'][0])*np.random.random())
    8
    9tc = pymc.Uniform('tc', bounds['tc'][0], bounds['tc'][1], \
    10     value = bounds['tc'][0] + (bounds['tc'][1]-bounds['tc'][0])*np.random.random())
    11
    12A = pymc.Uniform('A', bounds['A'][0], bounds['A'][1], \
    13    value = bounds['A'][0] + (bounds['A'][1]-bounds['A'][0])*np.random.random())
    14
    15dA = pymc.Uniform('dA', bounds['dA'][0], bounds['dA'][1], \
    16     value = bounds['dA'][0] + (bounds['dA'][1]-bounds['dA'][0])*np.random.random())
    17
    18f = pymc.Uniform('f', bounds['f'][0], bounds['f'][1], \
    19    value = bounds['f'][0] +  (bounds['f'][1]-bounds['f'][0])*np.random.random())
    20
    21df = pymc.Uniform('df', bounds['df'][0], bounds['df'][1], \
    22     value = bounds['df'][0] +  (bounds['df'][1]-bounds['df'][0])*np.random.random())

    然后定义我们的模板。在pymc中,这只是另一个变量,但是我们需要告诉它这不是一个随机变量,它的值取决于所有其它变量。我们使用deterministic装饰器来进行此操作。

    1@pymc.deterministic
    2def y_model(t=t, tc=tc, offset=offset, A=A, dA=dA, f=f, df=df):
    3    return chirp(t, tc, offset, A, dA, f, df)

    其次告诉pymc我们的数据。数据由pymc也表示为变量,这在一开始可能有点反直觉。我们使用pymc.Normal()来处理高斯噪声,该分布的均值为模板y_model,其方差由噪声强度给出。

    要指定它是数据而不是变量,我们使用两个额外的关键字来指定它:我们设置observed=True表示它是一个在采样过程中永不改变的观测值,我们使用value=y_obs将它的值设置为我们的实际数据。这使得y是一个随机变量实体,而不是一个要采样的变量。接下来,我们将所有内容包装到一个pymc模型中,并启动一个新的MCMC实例。

    1y = pymc.Normal('y', mu=y_model, tau=sigma**-2, observed=True, value=y_obs)
    2M = pymc.MCMC([y, tc, offset, A, dA, f, df, y_model])

    最后告诉pymc我们需要的样本数,它就会为我们做所有的事情。

    burn参数指定预烧期(burn-in period),即过多少步后马尔科夫链是趋向稳定分布的(通常是经验性的)。MCMC采样将花费一些时间来找到要采样参数空间的正确部分,而运行开始时的采样通常并不有趣,并且会扭曲我们对概率分布的描述,因此我们可以丢弃其中一些。

    thin参数指定如何稀释样本。MCMC样本有一些自相关,所以你可以只取第n个样本让它们相互独立。我们暂时不用担心,并将其设置为默认值1。

    1M.sample(iter=100000, burn=50000, thin=1)

    使用trace可以输出模型中每个变量的样本,例如M.trace('A')[:]将给出A的样本。

    为方便起见,我们将保存结果,然后将其加载回去,以便我们可以使用预先生成的数据。

     1# Convert the results of our MCMC sampling
    2def mcmc_dataframe(M, skip_vars):
    3    names = []
    4    data = []
    5    for var in M.variables:
    6        try:
    7            name = var.trace.name
    8            if name not in skip_vars:
    9                names.append(var.trace.name)
    10                data.append(var.trace.gettrace().tolist())
    11        except AttributeError:
    12            pass
    13    df = pd.DataFrame(data).T
    14    df.columns = names
    15    return df
    16
    17# Save as csv
    18mcmc_dataframe(M, ['y_model']).to_csv('samples.csv')

    03

    结果

    加载数据,绘制散点图矩阵(图中红叉为真实参数值)。

     1# Some settings for Seaborn
    2sns.set(context='notebook', style="ticks", color_codes=True)
    3
    4# Wrapper function for the hexbin plotting style
    5def hexbin(x, y, color, **kwargs):
    6    cmap = sns.light_palette(color, as_cmap=True)
    7    plt.hexbin(x, y, gridsize=15, cmap = cmap, mincnt = 1)
    8    return
    9
    10# Helper arrays
    11var_order = ['offset''tc''A''dA''f''df']
    12truth = [offset_true, tc_true, A_true, dA_true, f_true, df_true]
    13samples = pd.read_csv('samples.csv')
    14est = samples[var_order].max().values
    15
    16# Plot a pairwise grid of scatter plots and histograms
    17with sns.axes_style("white"):
    18    g = sns.pairplot(
    19        samples[var_order],
    20        # We don't actually want the scatter plots, so we set the marker size to 0
    21        plot_kws={"s"0}, 
    22        diag_kws={'color''green''lw'0})
    23
    24    # Replace the scatter plots with hexbin plots
    25    g.map_offdiag(hexbin, color='green');
    26
    27# Plot the injected values in the histograms
    28for i in range(len(var_order)):
    29    for j in range(len(var_order)):
    30        if i != j:
    31            ax = g.axes[i,j]
    32            ax.scatter([truth[j]], [truth[i]], color = "r", marker = 'x', s=100, linewidths=3)

    0d4db085439b30649fe0a3d4596b4f83.png

    上面的散点图矩阵让我们对结果能有一个很好的了解,它向我们展示了参数的估计情况,并清楚地说明了参数之间的相关性。

    我们看到对AdA的估算并不是很好,这不是我们技术的问题,而是数据根本不包含我们准确估计它们所需的信息,并且它们之间还存在负相关性。同样,fdf也显示出很强的负相关性,通常很难同时估计它们。然而,它们的每个峰值都非常接近注入值。并合时间tc的估计非常准确,在正确的值上出现了一个尖峰,其估计值与其它参数似乎都没有相关性。

    04

    结语

    引力波探测器面临的巨大挑战是信号要比噪声小几个数量级,但是由于可以根据广义相对论预测信号的形状,因此我们仍可以将信号从噪声中挖掘出来。

    参数估计的实际情况比我们这里看过的例子要复杂得多,它们通常包含更多的噪声和更多的参数,这意味着要进行更高维度的采样,需要更高级的技术,例如嵌套采样。但是,我希望这个笔记至少给出了如何使用匹配滤波来探测引力波的想法。

    点击文末左下角「

    来源:https://marcotompitak.github.io/pe-pymc/

    6c458169b10545aacb8f96f344e4c0a3.png

    附录:MCMC方法的代表算法--Metropolis-Hastings算法

    输入任意选定的建议分布(状态转移核)  ,抽样的目标分布密度数  ,收敛步数  要样本数  

    Step 1:随机选择一个初始值  ,样本集合 samples=[]

    Step 2:for t=0 to m+n:

    • 按照建议分布  随机抽取一个候选状态  

    • 计算接受概率: 

    • 从区间  中按均匀分布随机抽取一个数  ,若  ,则状态  ,否则,  保持不变

    • if (t >= m) 将  加入到 samples 中

    Step 3返回样本集合 samples。

    一个简单的Python例子如下,可以看出用MCMC采样得到的分布与真实分布相一致。

      1import random
    2import numpy as np
    3import matplotlib.pyplot as plt
    4
    5def mh(q, p, m, n):
    6    # randomize a number
    7    x = random.uniform(0.11)
    8    for t in range(0, m+n):
    9        x_sample = q.sample(x)
    10        try:
    11            accept_prob = min(1, p.prob(x_sample)*q.prob(x_sample, x)/(p.prob(x)*q.prob(x, x_sample)))
    12        except:
    13            accept_prob = 0
    14
    15        u = random.uniform(01)
    16
    17        if u  18            x = x_sample
    19
    20        if t >= m:
    21            yield x
    22
    23class Exponential(object):
    24    def __init__(self, scale):
    25        self.scale = scale
    26        self.lam = 1.0 / scale
    27
    28    def prob(self, x):
    29        if x <= 0:
    30            raise Exception("The sample shouldn't be less than zero")
    31
    32        result = self.lam * np.exp(-x * self.lam)
    33        return result
    34
    35    def sample(self, num):
    36        sample = np.random.exponential(self.scale, num)
    37        return sample
    38
    39# 假设我们的目标概率密度函数p1(x)是指数概率密度函数
    40scale = 5
    41p1 = Exponential(scale)
    42
    43
    44class Norm():
    45    def __init__(self, mean, std):
    46        self.mean = mean
    47        self.sigma = std
    48
    49    def prob(self, x):
    50        return np.exp(-(x - self.mean) ** 2 / (2 * self.sigma ** 2.0)) * 1.0 / (np.sqrt(2 * np.pi) * self.sigma)
    51
    52    def sample(self, num):
    53        sample = np.random.normal(self.mean, self.sigma, size=num)
    54        return sample
    55
    56# 假设我们的目标概率密度函数p1(x)是均值标准差分别为3,2的正态分布
    57p2 = Norm(32)
    58
    59class Transition():
    60    def __init__(self, sigma):
    61        self.sigma = sigma
    62
    63    def sample(self, cur_mean):
    64        cur_sample = np.random.normal(cur_mean, scale=self.sigma, size=1)[0]
    65        return cur_sample
    66
    67    def prob(self, mean, x):
    68        return np.exp(-(x-mean)**2/(2*self.sigma**2.0)) * 1.0/(np.sqrt(2 * np.pi)*self.sigma)
    69
    70# 假设我们的转移核标准差为10的正态分布
    71q = Transition(10)
    72
    73m = 100
    74n = 100000 # 采样个数
    75
    76simulate_samples_p1 = [li for li in mh(q, p1, m, n)]
    77
    78plt.figure(figsize=(88)) 
    79
    80plt.subplot(2,2,1)
    81plt.hist(simulate_samples_p1, 100, color='blue')
    82plt.title("Simulated X ~ Exponential(1/5)")
    83
    84samples = p1.sample(n)
    85plt.subplot(2,2,2)
    86plt.hist(samples, 100, color='yellowgreen')
    87plt.title("True X ~ Exponential(1/5)")
    88
    89simulate_samples_p2 = [li for li in mh(q, p2, m, n)]
    90plt.subplot(2,2,3)
    91plt.hist(simulate_samples_p2, 50, color='blue')
    92plt.title("Simulated X ~ N(3,2)")
    93
    94samples = p2.sample(n)
    95plt.subplot(2,2,4)
    96plt.hist(samples, 50, color='yellowgreen')
    97plt.title("True X ~ N(3,2)")
    98
    99plt.suptitle("Transition Kernel N(0,10) simulation results")
    100plt.show()

    c00b7087fdb0ec5d1ed027e960f0ce45.png

    来源:https://zhuanlan.zhihu.com/p/108621943

    70dbec6b863ccf1c87a09e0cf2e483fa.png

    长按下方图片关注「引力天文」查看更多历史文章

    b53daeb58ffec9ca2467f467e40963d8.png

    70dbec6b863ccf1c87a09e0cf2e483fa.png留言请直接在公众号里发消息
    展开全文
  • 参数估计的计算方法

    千次阅读 2020-05-27 19:21:58
    参数估计的计算方法极大后验(MAP)及拉普拉斯逼近基于马尔可夫链的蒙特卡洛参数推断(MCMC)期望极大化(EM) (参数估计所有内容) 极大后验(MAP)及拉普拉斯逼近 极大后验估计: MAP是通过确定后验分布的极大值得到的,...

    参数估计所有内容

    极大后验(MAP)及拉普拉斯逼近

    极大后验估计:
    MAP是通过确定后验分布的极大值得到的,在点估计中的表达式为:在这里插入图片描述MAP 估计可等效为能量函数的极小值:
    在这里插入图片描述其中,能量函数表达式为:
    在这里插入图片描述参数的极大似然估计即为一种标准先验概率满足下式的MAP估计:
    在这里插入图片描述

    能量函数的极小值可通过多种基于无梯度或梯度的广义最优算法计算得到。
    若要利用基于梯度的广义最优化算法,需要计算能量函数的导数。计算能量函数导数有两种方法:
    1). 对能量函数的递归方程用一种特殊的方法进行微分,所得的结果称为灵敏度函数,可由类似于滤波过程中递归方程计算得到。
    2). 利用Fisher特性,将能量函数的梯度表示为平滑分布所有数据对数似然导数的期望值。该方法相比于直接微分法优点在于不要计算额外的递归方程。

    MAP估计的缺点:该方法实际运用了狄拉克 δ\bm{\delta} 方程逼近后验分布,即:
    在这里插入图片描述从而忽略了后验分布的变化。

    拉普拉斯逼近法:
    利用能量函数的二阶导数(Hessian) 实现对后验分布的高斯逼近,表达式为:
    在这里插入图片描述式中,H(θ^MAP)\bm{H(\hat\theta^{MAP})}MAP 估计值出的 Hessian 矩阵。该方法需要计算或者逼近能量函数的二阶导数。

    基于马尔可夫链的蒙特卡洛参数推断(MCMC)

    MH算法:
    Metropolis-Hastings(MH)算法是最为常见的一种MCMC方法。MH引入建议密度 q(θ(i)θ(i1))\bm{q(\theta^{(i)}|\theta^{(i-1)})},通过结合已知采样点 θ(i1)\bm{\theta^{(i-1)}},给出新的建议采样点 θ(i)\bm{\theta^{(i)}}
    MH算法步骤为:
    首选由任意初始分布得到起始点 θ(0)\bm{\theta^{(0)}}
    然后,对于 i=1,2,,N\bm{i=1,2,\cdot\cdot\cdot,N},有

    1). 从建议密度中选择一个参考点 θ\bm{\theta^{*}}:
    在这里插入图片描述2). 计算接受概率:
    在这里插入图片描述3). 产生均匀分布随机变量 uU(0,1)\bm{u \sim U(0,1)},并令:
    在这里插入图片描述

    Metropolis算法为Metropolis-Hastings算法的一种特例,其建议分布为对称的,即 q(θ(i)θ(i1))=q(θ(i1)θ(i))\bm{q(\theta^{(i)}|\theta^{(i-1)})}=\bm{q(\theta^{(i-1)}|\theta^{(i)})},此时,接受概率缩减为:
    在这里插入图片描述
    建议分布的选取对Metropolis-Hastings算法的性能影响较大,通常采用高斯部分作为建议分布。即:
    在这里插入图片描述式中,i1\bm{\sum^{i-1}} 为合适的协方差矩阵。
    选定协方差矩阵其中一种方法是自适应马尔可夫链蒙特卡洛法(AMCMC\bm{AMCMC}),该方法可以载解算MCMC\bm{MCMC}时,自动调整所需要的高斯建议分布协方差阵。
    AMCMC\bm{AMCMC} 的思想为:利用先前采样点产生的协方差阵作为后验分布实际协方差阵的估计值。在已知先前的协方差矩阵的情况下,就可以计算建议分布的协方差阵。

    RAM(robust adaptive Metropolis)算法:
    RAM步骤为:
    1). 有初始分布 p0(θ)\bm{p_0(\theta)} 获得 θ(0)\bm{\theta^{(0)}},利用初始协方差阵的下三角阵 Chollesky\bm{Chollesky} 因子对 S0\bm{S_0} 进行初始化。
    2). 获取一个采样点 θ=θi1+Si1ri,riN(0,I)\bm{\theta^*=\theta_{i-1}+S_{i-1}r_i},r_i\sim N(0,I)
    3). 计算接受概率:
    在这里插入图片描述4). 从均匀分布 U(0,1)\bm{U(0,1)} 中采样的到一个服从均匀分布的随机变量 u\bm{u}
    5).
    在这里插入图片描述6). 计算下三角矩阵 Si\bm{S_i},使得该矩阵满足对角线元素为正,且该阵满足:
    在这里插入图片描述式中,{η}i1(0,1]\bm{\{\eta\}_{i\ge1}\subset(0,1]},为一个自适应步长序列,并逐渐衰减为零。可任意选取,在其他文献中给出一个选取建议:ηi=iγ,γ(1/2,1]\bm{\eta_i=i^{- \gamma},\gamma\in(1/2,1]}
    7).ii+1\bm{i\leftarrow i+1},跳至步骤2),直到获取足够多的采样点。

    期望极大化(EM)

    EM\bm{EM} 算法适用于边缘似然无法计算但仍能计算得到似然函数下界的情形。设 q(x0:T)\bm{q(x_{0:T})} 为状态的任意概率密度函数,则有:
    在这里插入图片描述式中,函数 F\bm{F} 的定义为:
    在这里插入图片描述
    EM\bm{EM} 算法中国最重要的就是通过迭代最大化下界 F[q(x0:T),θ]\bm{F[q(x_{0:T}),\theta]},进而实现 logp(y1:Tθ)\bm{logp(y_{1:T}|\theta)} 的最大化。

    简易的EM算法
    对函数 F\bm{F} 的最大化,可通过协调提升下述步骤实现:
    假设初始值为:q(0),θ(0)\bm{q^{(0)},\theta^{(0)}}
    对于 n=0,1,2,\bm{n=0,1,2,\cdot\cdot\cdot},执行如下步骤:
    1). E-步骤:获取:
    在这里插入图片描述2). M-步骤:获取:
    在这里插入图片描述
    为实现 EM\bm{EM} 算法,上式的极大化必行是可行的,所幸在其他文献中给出 E\bm{E} 步骤极大化后的结果:
    在这里插入图片描述将上式带入 F\bm{F} 函数中,得:
    在这里插入图片描述
    由于在实行 M\bm{M} 步骤中,上式右边第二项与 θ\bm{\theta} 无关,故极大化 函数 F\bm{F} 就得等价于对等式右边第一项对极大化,因此,在 EM\bm{EM} 算法中,可表示为:
    在这里插入图片描述上式为 EM\bm{EM} 算法得到关于 p(x0:T,y0:Tθ)\bm{p(x_{0:T},y_{0:T}|\theta)} 的期望,p(x0:T,y0:Tθ)\bm{p(x_{0:T},y_{0:T}|\theta)}是在参数 θ\bm{\theta} 条件下关于全部状态量和量测量的联合后验似然分布。

    EM\bm{EM} 算法:
    假设初始值为 θ(0)\bm{\theta^{(0)}}。对于 n=0,1,2,\bm{n=0,1,2,\cdot\cdot\cdot},执行如下步骤:
    1). E-步骤:计算 Q(θ,θ(n))\bm{\mathcal{Q(\theta,\theta^{(n)})}}
    2). M-步骤:计算:
    在这里插入图片描述结合下式:状态空间模型的马尔可夫链结构:
    在这里插入图片描述完全数据的对数似然表达式为:
    在这里插入图片描述因此,Q\bm{\mathcal{Q}} 的表达式及上述算法的 E\bm{E} 步骤可以简化为:
    在这里插入图片描述式中:
    在这里插入图片描述上式均为关于平滑分布的期望,即可不必在计算完全后验分布的期望。

    EM\bm{EM} 算法的 E\bm{E} 步骤中,需要对 Q\bm{\mathcal{Q}} 的表达式关于参数 θ\bm{\theta} 的极大化。
    实用的方法为将梯度设为零,求的此时极大值时候的 θ\bm{\theta} 取值:
    在这里插入图片描述
    EM\bm{EM} 算法,
    E\bm{E} 步骤:利用上一步得到的表达式 Q(θ,θn)\bm{\mathcal{Q(\theta,\theta^n)}} 求出满足表达式 Q(θ,θn)\bm{\mathcal{Q(\theta,\theta^n)}} 极大值的 参数 θn\bm{\theta^n}
    M\bm{M} 步骤:将 E\bm{E} 步骤求出的参数 θn\bm{\theta^n} 作为下一步已知的值 θn+1\bm{\theta^{n+1}} 带入表达式 Q(θ,θn)\bm{\mathcal{Q(\theta,\theta^n)}} 中得到Q(θ,θn+1)\bm{\mathcal{Q(\theta,\theta^{n+1})}},进而最大化函数 F\bm{F} 的下限。
    循环重复 E,M\bm{E,M} 步骤,直到满足要求。

    展开全文
  • 文章目录Zero.写作动机一、模型原理二、编程实现Notice. ...网络上关于蒙特卡洛方法几乎清一色都是在介绍Buffon实验并以此估计某个量。这里,我们介绍蒙特卡洛树用于参数寻优。 一、模型原理 下面推荐几个博客,这些...

    Zero.写作动机

    对给定参数区间内部进行搜索,寻找到最优参数近似解的方法有很多。比如网格搜索。但是网格搜索太过暴力,往往花销过大。这里介绍一种新的参数寻优方法——蒙特卡洛树搜索
    网络上关于蒙特卡洛方法几乎清一色都是在介绍Buffon实验并以此估计某个量。这里,我们介绍蒙特卡洛树用于参数寻优。

    一、模型原理

    下面推荐几个博客,这些文章已经介绍得很好了:

    ①https://blog.csdn.net/ljyt2/article/details/78332802
    ②https://www.jianshu.com/p/a34f06885ef8

    二、编程实现

    Version one: Python
    https://www.jianshu.com/p/a34f06885ef8

    Version Two: Matlab
    鉴于实际需求,笔者在Python版本的基础上实现了matlab版本,涉及到matlab的面向对象编程。读者诸君按需获取即可

    state.m文件

    classdef State < handle
        properties
            value
            round
            choices
            PATH
            x2 %因为假定现在只用MCST找到第三步迭代的最优参数
            y 
            sigma
            im
        end
        methods
            function self= State(x2,y,sigma,im)
            %在这里进行初始化
            self.value = 0;
            self.round = 0;
            self.choices = [];
            self.PATH = [0.1:0.2:3];
            self.x2 = x2;
            self.y = y;
            self.sigma = sigma;
            self.im = im;
           
            end
            
            function state = new_state(self)
                choice = randperm(numel(self.PATH));
                choice = self.PATH(choice(1));%从一维数组中进行随机采样
                state = State(self.x2,self.y,self.sigma,self.im);
                %对于辣椒的彩色图片,第三步迭代的默认两个参数是0.7, 0.8
                value_ = 0;
                if numel(self.choices) == 1 %当前在选择第二个参数
                    %计算潜在的value
                    x3 = step(self.x2, self.y, self.sigma^2, 15, 7, self.choices(1), choice);
                    value_ = - (sum(sum((x3 - self.im).^2)) / numel(x3)); %反向来
                elseif numel(self.choices) == 0 %当前在选择第一个参数
                    %计算潜在的value
                    x3 = step(self.x2, self.y, self.sigma^2, 15, 7, choice,0.8);
                    value_ = - (sum(sum((x3 - self.im).^2)) / numel(x3)); %反向来
                else
                    value_ = 0;
                end
                
                
                %得到一个参数的选择结果
                state.value = self.value +  value_; %价值计算函数需要更改
                state.round = self.round+1;
                state.choices = [ self.choices,choice ];%扩充当前的选择
            end
            
            function display(self)
                fprintf(1,'class State:\n');%表示在终端上进行输出
                fprintf(1,'value = %f\n',self.value);
                fprintf(1,'round = %d\n',self.round);
                fprintf(1,'ready to show the choice array:\n');
                for i = 1:numel(self.choices)
                    if i == 1
                        fprintf(1,'[');
                    end
                    fprintf(1,'%d,',self.choices(i));
                    if i == numel(self.choices)
                        fprintf(1,']');
                    end
                end
            end
            
        end
    end
    

    Node.m文件

    classdef Node < handle
        properties
            parent
            children 
            quality
            visit
            state
            MAX_DEPTH = 2
            MAX_CHOICE = numel([0.1:0.2:3]) %其实代表的是children数组的长度的上限
        end
        methods
            function self= Node()
                self.quality = 0.0;
                self.visit = 0;
               %剩下的变量没有定义
            end
            
            function add_child(self,node)
                
                fprintf(1,'printing node in function add_child\n');
                node          
                
                self.children = [self.children,node];
                node.parent = self;
            end
            
            function display(self)
                fprintf(1,'class Node:\n');%表示在终端上进行输出
                fprintf(1,'quality = %f\n',self.quality);
                fprintf(1,'visit = %d\n',self.visit);
            end
            
            function  child_node = expand(cnt_node)
                %随机选择一个之前没有扩展过的——也就是不在children列表中的一个子节点进行扩展,随机性在new_state的时候的随机函数中体现出来
                %返回当前结点扩展出的子节点
                fprintf(1,'printing node in function EXPAND\n');
                cnt_node
                fprintf(1,'printing value of the ori_state in function EXPAND:%f\n\n',cnt_node.state.value);
                cnt_node.state.choices
                
                state = new_state(cnt_node.state);
                %拿到当前结点的children列表中的子节点的状态
                sub_state_value_list = [];
                for i = 1:numel(cnt_node.children)
                    sub_state_value_list(i) = cnt_node.children(i).state.value;
                end
                fprintf(1,'printing value of the new_state in function EXPAND:%f\n\n',state.value);
                state.choices
                
                while ismember(state.value,sub_state_value_list)
                    fprintf(1,'printing value of the new_state in function EXPAND:%f\n\n',state.value);
                     state.choices
                    state = new_state(cnt_node.state);
                end
                child_node = Node();
                child_node.state = state;
                add_child(cnt_node,child_node);
                
                
                fprintf(1,'printing value of the end_child_state in function EXPAND\n');            
                for i = 1:numel(cnt_node.children)
                    fprintf(1,'printing value of the end_child_state in function EXPAND:%f\n\n',cnt_node.children(i).state.value);
                    cnt_node.children(i).state.choices
                end                        
                
            end
            
            function best = best_child(node)
                %返回当前结点的children列表中最适合作为扩展结点的子节点
                
                fprintf(1,'printing node in function BEST_CHILD\n');
                node   
                
                best_score = -100000000; %代表负无穷
                best = -1 ;%初始化
                for i=  1:numel(node.children)
                    C = 1/sqrt(2.0);
                    sub_node = node.children(i);
                    left = sub_node.quality / sub_node.visit; %分母是被访问的次数
                    right = 2.0*log(node.visit)/sub_node.visit;
                    score = left+C*sqrt(right);
                    
                    if score >best_score
                        best = sub_node;
                        best_score = score;
                    end
                end
            end
            
            function node = tree_policy(node)
                fprintf(1,'printing node in function TREE_POLICY\n');
                node   
            %选择+expand扩展
            %调用逻辑:如果当前结点还有子节点没有被添加到children列表——也就是还没有expand过,那么就从还没有扩展过的子节点中随机选择一个进行扩展,并返回该被需选中的子节点
            %调用逻辑:如果当前结点是叶子结点,直接返回该结点
            %调用逻辑:如果当前结点的所有子节点都已经被加入到了children列表,那么就从中选择一个收益最高的结点进行扩展,并且返回该结点
                %选择是否是叶子结点
                count = 0;
                while node.state.round < node.MAX_DEPTH
                    fprintf(1,'running while-end with count:%d in Node.m/line73\n',count);
                    count = count +1;
                    if numel(node.children) < node.MAX_CHOICE
                        node = expand(node);
                        return
                    else
                        node = best_child(node);
                    end
                end     
            end
            
            function expanded_value = default_policy(node)
                fprintf(1,'printing node in function DEFAULT_POLICY\n');
                node   
            %模拟
            %算一次从当前结点随机走到叶节点的收益
                now_state = node.state;
                count= 0;
                while now_state.round < node.MAX_DEPTH
                    fprintf(1,'running while-end with count:%d in Node.m/line90\n',count);
                    count = count +1;
                    now_state = new_state(now_state);
                end
                
                expanded_value = now_state.value;
                
            end
            
            function backup(node,reward)
                fprintf(1,'printing node in function BACKUP\n');
                node   
                %从当前结点带着reward回溯到根节点,并且增加路径上的每个结点的visit次数和quality
                while ~isempty(node)               
                    fprintf(1,'not empty\n');
                    node.visit = node.visit +1;
                    node.quality = node.quality+reward;
                    node = node.parent;
                end
            
            end
            
            function best = mcts(node)
                %似乎是多次尝试扩展,选择当前扩展到children列表中的子节点中的收益最好的一个子结点进行扩展,并且返回该被选中的子节点
              %  times =  5 ;%为什么是5?
                times = 20;
                for i = 1:times
                    expand = tree_policy(node);%当前结点向下选择扩展一个结点
                    reward = default_policy(expand);%计算从该扩展结点走到叶子结点的随机一条路径的一种收益情况
                    backup(expand,reward);
                end
                best = best_child(node);
                
            end
            
            function main(self)
                init_state = State();
                init_node = Node();
                init_node.state = init_state;
                cnt_node = init_node;
                
                
                for i = 1:self.MAX_DEPTH
                    cnt_node = mcts(cnt_node);
                end
                
            end
            
        end
    end
    
    

    Notice.

    在matlab的实现版本中,注意两种不同的类的写法classdef name < handle是引用类型,这样的类可以作为另外一个类的属性存在。classdef name是按value类型,这样的类如果想要使用自己的实例对象作为类的一个属性会报错。
    上面的Node类和State类都属于引用类型。

    展开全文
  • 对含有MA(3)误差项的面板数据模型,应用半参数计量经济学理论,在相关假设基础上推导出半参数估计方法,并根据蒙特卡洛模拟对这种方法进行实证检验,结果证明这种半参数估计方法是可行的,并且相对参数估计方法(最小二乘...
  • 蒙特卡洛模拟法简介蒙特卡洛(MonteCarlo)模拟是一种通过设定随机过程,反复生成时间序列,计算参数估计量和统计量,进而研究其分布特征的方法。具体的,当系统中各个单元的可靠性特征量已知,但系统的可靠性过于...

    一 蒙特卡洛模拟法简介

    蒙特卡洛(Monte

    Carlo)模拟是一种通过设定随机过程,反复生成时间序列,计算参数估计量和统计量,进而研究其分布特征的方法。具体的,当系统中各个单元的可靠性特征量已知,但系统的可靠性过于复杂,难以建立可靠性预计的精确数学模型或模型太复杂而不便应用时,可用随机模拟法近似计算出系统可靠性的预计值;随着模拟次数的增多,其预计精度也逐渐增高。由于涉及到时间序列的反复生成,蒙特卡洛模拟法是以高容量和高速度的计算机为前提条件的,因此只是在近些年才得到广泛推广。

    这个术语是二战时期美国物理学家Metropolis执行曼哈顿计划的过程中提出来的。

    蒙特卡洛模拟方法的原理是当问题或对象本身具有概率特征时,可以用计算机模拟的方法产生抽样结果,根据抽样计算统计量或者参数的值;随着模拟次数的增多,可以通过对各次统计量或参数的估计值求平均的方法得到稳定结论。

    二 蒙特卡洛模拟法求解步骤

    应用此方法求解工程技术问题可以分为两类:确定性问题和随机性问题。

    解题步骤如下:

    1.根据提出的问题构造一个简单、适用的概率模型或随机模型,使问题的解对应于该模型中随机变量的某些特征(如概率、均值和方差等),所构造的模型在主要特征参量方面要与实际问题或系统相一致

    2

    .根据模型中各个随机变量的分布,在计算机上产生随机数,实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数,然后生成服从某一分布的随机数,方可进行随机模拟试验。

    3.

    根据概率模型的特点和随机变量的分布特性,设计和选取合适的抽样方法,并对每个随机变量进行抽样(包括直接抽样、分层抽样、相关抽样、重要抽样等)。

    4.按照所建立的模型进行仿真试验、计算,求出问题的随机解。

    5. 统计分析模拟试验结果,给出问题的概率解以及解的精度估计。

    三 蒙特卡洛模拟法的应用领域

    蒙特卡洛模拟法的应用领域主要有:

    1.直接应用蒙特卡洛模拟:应用大规模的随机数列来模拟复杂系统,得到某些参数或重要指标。

    2.蒙特卡洛积分:利用随机数列计算积分,维数越高,积分效率越高。

    3.MCMC:这是直接应用蒙特卡洛模拟方法的推广,该方法中随机数的产生是采用的马尔科夫链形式。

    四 资产组合模拟

    假设有五种资产,其日收益率(%)分别为

    0.0246 0.0189 0.0273 0.0141 0.0311

    标准差分别为

    0.9509 1.4259, 1.5227, 1.1062, 1.0877

    相关系数矩阵为

    1.0000 0.4403 0.4735 0.4334 0.6855

    0.4403 1.0000 0.7597 0.7809 0.4343

    0.4735 0.7597 1.0000 0.6978 0.4926

    0.4334 0.7809 0.6978 1.0000 0.4289

    0.6855 0.4343 0.4926 0.4289 1.0000

    假设初始价格都为100,模拟天数为504天,模拟线程为2,程序如下

    %run.m

    ExpReturn = [0.0246 0.0189 0.0273 0.0141 0.0311]/100; %期望收益

    Sigmas = [0.9509 1.4259, 1.5227, 1.1062, 1.0877]/100;%标准差

    Correlations = [1.0000 0.4403 0.4735 0.4334 0.6855

    0.4403 1.0000 0.7597 0.7809 0.4343

    0.4735 0.7597 1.0000 0.6978 0.4926

    0.4334 0.7809 0.6978 1.0000 0.4289

    0.6855 0.4343 0.4926 0.4289 1.0000

    ];%相关系数

    ExpCov = corr2cov(Sigmas, Correlations);%协方差

    StartPrice = 100;%初始价格

    NumObs = 504;

    NumSim = 2;

    RetIntervals = 1;

    NumAssets = 5;

    %开始模拟

    randn('state', 0);

    RetExact = portsim(ExpReturn, ExpCov, NumObs, RetIntervals,

    NumSim);

    Weights = ones(NumAssets, 1)/ NumAssets;

    PortRetExact = zeros(NumObs, NumSim);

    for i = 1:NumSim

    PortRetExact(:, i) = RetExact(:,:,i)*Weights;

    end

    PortExact = ret2tick(PortRetExact,

    repmat(StartPrice, 1, NumSim));

    plot(PortExact, '-r');

    一 蒙特卡洛模拟法简介

    蒙特卡洛(Monte

    Carlo)模拟是一种通过设定随机过程,反复生成时间序列,计算参数估计量和统计量,进而研究其分布特征的方法。具体的,当系统中各个单元的可靠性特征量已知,但系统的可靠性过于复杂,难以建立可靠性预计的精确数学模型或模型太复杂而不便应用时,可用随机模拟法近似计算出系统可靠性的预计值;随着模拟次数的增多,其预计精度也逐渐增高。由于涉及到时间序列的反复生成,蒙特卡洛模拟法是以高容量和高速度的计算机为前提条件的,因此只是在近些年才得到广泛推广。

    这个术语是二战时期美国物理学家Metropolis执行曼哈顿计划的过程中提出来的。

    蒙特卡洛模拟方法的原理是当问题或对象本身具有概率特征时,可以用计算机模拟的方法产生抽样结果,根据抽样计算统计量或者参数的值;随着模拟次数的增多,可以通过对各次统计量或参数的估计值求平均的方法得到稳定结论。

    二 蒙特卡洛模拟法求解步骤

    应用此方法求解工程技术问题可以分为两类:确定性问题和随机性问题。

    解题步骤如下:

    1.根据提出的问题构造一个简单、适用的概率模型或随机模型,使问题的解对应于该模型中随机变量的某些特征(如概率、均值和方差等),所构造的模型在主要特征参量方面要与实际问题或系统相一致

    2

    .根据模型中各个随机变量的分布,在计算机上产生随机数,实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数,然后生成服从某一分布的随机数,方可进行随机模拟试验。

    3.

    根据概率模型的特点和随机变量的分布特性,设计和选取合适的抽样方法,并对每个随机变量进行抽样(包括直接抽样、分层抽样、相关抽样、重要抽样等)。

    4.按照所建立的模型进行仿真试验、计算,求出问题的随机解。

    5. 统计分析模拟试验结果,给出问题的概率解以及解的精度估计。

    三 蒙特卡洛模拟法的应用领域

    蒙特卡洛模拟法的应用领域主要有:

    1.直接应用蒙特卡洛模拟:应用大规模的随机数列来模拟复杂系统,得到某些参数或重要指标。

    2.蒙特卡洛积分:利用随机数列计算积分,维数越高,积分效率越高。

    3.MCMC:这是直接应用蒙特卡洛模拟方法的推广,该方法中随机数的产生是采用的马尔科夫链形式。

    四 资产组合模拟

    假设有五种资产,其日收益率(%)分别为

    0.0246 0.0189 0.0273 0.0141 0.0311

    标准差分别为

    0.9509 1.4259, 1.5227, 1.1062, 1.0877

    相关系数矩阵为

    1.0000 0.4403 0.4735 0.4334 0.6855

    0.4403 1.0000 0.7597 0.7809 0.4343

    0.4735 0.7597 1.0000 0.6978 0.4926

    0.4334 0.7809 0.6978 1.0000 0.4289

    0.6855 0.4343 0.4926 0.4289 1.0000

    假设初始价格都为100,模拟天数为504天,模拟线程为2,程序如下

    %run.m

    ExpReturn = [0.0246 0.0189 0.0273 0.0141 0.0311]/100; %期望收益

    Sigmas = [0.9509 1.4259, 1.5227, 1.1062, 1.0877]/100;%标准差

    Correlations = [1.0000 0.4403 0.4735 0.4334 0.6855

    0.4403 1.0000 0.7597 0.7809 0.4343

    0.4735 0.7597 1.0000 0.6978 0.4926

    0.4334 0.7809 0.6978 1.0000 0.4289

    0.6855 0.4343 0.4926 0.4289 1.0000

    ];%相关系数

    ExpCov = corr2cov(Sigmas, Correlations);%协方差

    StartPrice = 100;%初始价格

    NumObs = 504;

    NumSim = 2;

    RetIntervals = 1;

    NumAssets = 5;

    %开始模拟

    randn('state', 0);

    RetExact = portsim(ExpReturn, ExpCov, NumObs, RetIntervals,

    NumSim);

    Weights = ones(NumAssets, 1)/ NumAssets;

    PortRetExact = zeros(NumObs, NumSim);

    for i = 1:NumSim

    PortRetExact(:, i) = RetExact(:,:,i)*Weights;

    end

    PortExact = ret2tick(PortRetExact,

    repmat(StartPrice, 1, NumSim));

    plot(PortExact, '-r');

    一 蒙特卡洛模拟法简介

    蒙特卡洛(Monte

    Carlo)模拟是一种通过设定随机过程,反复生成时间序列,计算参数估计量和统计量,进而研究其分布特征的方法。具体的,当系统中各个单元的可靠性特征量已知,但系统的可靠性过于复杂,难以建立可靠性预计的精确数学模型或模型太复杂而不便应用时,可用随机模拟法近似计算出系统可靠性的预计值;随着模拟次数的增多,其预计精度也逐渐增高。由于涉及到时间序列的反复生成,蒙特卡洛模拟法是以高容量和高速度的计算机为前提条件的,因此只是在近些年才得到广泛推广。

    这个术语是二战时期美国物理学家Metropolis执行曼哈顿计划的过程中提出来的。

    蒙特卡洛模拟方法的原理是当问题或对象本身具有概率特征时,可以用计算机模拟的方法产生抽样结果,根据抽样计算统计量或者参数的值;随着模拟次数的增多,可以通过对各次统计量或参数的估计值求平均的方法得到稳定结论。

    二 蒙特卡洛模拟法求解步骤

    应用此方法求解工程技术问题可以分为两类:确定性问题和随机性问题。

    解题步骤如下:

    1.根据提出的问题构造一个简单、适用的概率模型或随机模型,使问题的解对应于该模型中随机变量的某些特征(如概率、均值和方差等),所构造的模型在主要特征参量方面要与实际问题或系统相一致

    2

    .根据模型中各个随机变量的分布,在计算机上产生随机数,实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数,然后生成服从某一分布的随机数,方可进行随机模拟试验。

    3.

    根据概率模型的特点和随机变量的分布特性,设计和选取合适的抽样方法,并对每个随机变量进行抽样(包括直接抽样、分层抽样、相关抽样、重要抽样等)。

    4.按照所建立的模型进行仿真试验、计算,求出问题的随机解。

    5. 统计分析模拟试验结果,给出问题的概率解以及解的精度估计。

    三 蒙特卡洛模拟法的应用领域

    蒙特卡洛模拟法的应用领域主要有:

    1.直接应用蒙特卡洛模拟:应用大规模的随机数列来模拟复杂系统,得到某些参数或重要指标。

    2.蒙特卡洛积分:利用随机数列计算积分,维数越高,积分效率越高。

    3.MCMC:这是直接应用蒙特卡洛模拟方法的推广,该方法中随机数的产生是采用的马尔科夫链形式。

    四 资产组合模拟

    假设有五种资产,其日收益率(%)分别为

    0.0246 0.0189 0.0273 0.0141 0.0311

    标准差分别为

    0.9509 1.4259, 1.5227, 1.1062, 1.0877

    相关系数矩阵为

    1.0000 0.4403 0.4735 0.4334 0.6855

    0.4403 1.0000 0.7597 0.7809 0.4343

    0.4735 0.7597 1.0000 0.6978 0.4926

    0.4334 0.7809 0.6978 1.0000 0.4289

    0.6855 0.4343 0.4926 0.4289 1.0000

    假设初始价格都为100,模拟天数为504天,模拟线程为2,程序如下

    %run.m

    ExpReturn = [0.0246 0.0189 0.0273 0.0141 0.0311]/100; %期望收益

    Sigmas = [0.9509 1.4259, 1.5227, 1.1062, 1.0877]/100;%标准差

    Correlations = [1.0000 0.4403 0.4735 0.4334 0.6855

    0.4403 1.0000 0.7597 0.7809 0.4343

    0.4735 0.7597 1.0000 0.6978 0.4926

    0.4334 0.7809 0.6978 1.0000 0.4289

    0.6855 0.4343 0.4926 0.4289 1.0000

    ];%相关系数

    ExpCov = corr2cov(Sigmas, Correlations);%协方差

    StartPrice = 100;%初始价格

    NumObs = 504;

    NumSim = 2;

    RetIntervals = 1;

    NumAssets = 5;

    %开始模拟

    randn('state', 0);

    RetExact = portsim(ExpReturn, ExpCov, NumObs, RetIntervals,

    NumSim);

    Weights = ones(NumAssets, 1)/ NumAssets;

    PortRetExact = zeros(NumObs, NumSim);

    for i = 1:NumSim

    PortRetExact(:, i) = RetExact(:,:,i)*Weights;

    end

    PortExact = ret2tick(PortRetExact,

    repmat(StartPrice, 1, NumSim));

    plot(PortExact, '-r');

    一 蒙特卡洛模拟法简介

    蒙特卡洛(Monte

    Carlo)模拟是一种通过设定随机过程,反复生成时间序列,计算参数估计量和统计量,进而研究其分布特征的方法。具体的,当系统中各个单元的可靠性特征量已知,但系统的可靠性过于复杂,难以建立可靠性预计的精确数学模型或模型太复杂而不便应用时,可用随机模拟法近似计算出系统可靠性的预计值;随着模拟次数的增多,其预计精度也逐渐增高。由于涉及到时间序列的反复生成,蒙特卡洛模拟法是以高容量和高速度的计算机为前提条件的,因此只是在近些年才得到广泛推广。

    这个术语是二战时期美国物理学家Metropolis执行曼哈顿计划的过程中提出来的。

    蒙特卡洛模拟方法的原理是当问题或对象本身具有概率特征时,可以用计算机模拟的方法产生抽样结果,根据抽样计算统计量或者参数的值;随着模拟次数的增多,可以通过对各次统计量或参数的估计值求平均的方法得到稳定结论。

    二 蒙特卡洛模拟法求解步骤

    应用此方法求解工程技术问题可以分为两类:确定性问题和随机性问题。

    解题步骤如下:

    1.根据提出的问题构造一个简单、适用的概率模型或随机模型,使问题的解对应于该模型中随机变量的某些特征(如概率、均值和方差等),所构造的模型在主要特征参量方面要与实际问题或系统相一致

    2

    .根据模型中各个随机变量的分布,在计算机上产生随机数,实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数,然后生成服从某一分布的随机数,方可进行随机模拟试验。

    3.

    根据概率模型的特点和随机变量的分布特性,设计和选取合适的抽样方法,并对每个随机变量进行抽样(包括直接抽样、分层抽样、相关抽样、重要抽样等)。

    4.按照所建立的模型进行仿真试验、计算,求出问题的随机解。

    5. 统计分析模拟试验结果,给出问题的概率解以及解的精度估计。

    三 蒙特卡洛模拟法的应用领域

    蒙特卡洛模拟法的应用领域主要有:

    1.直接应用蒙特卡洛模拟:应用大规模的随机数列来模拟复杂系统,得到某些参数或重要指标。

    2.蒙特卡洛积分:利用随机数列计算积分,维数越高,积分效率越高。

    3.MCMC:这是直接应用蒙特卡洛模拟方法的推广,该方法中随机数的产生是采用的马尔科夫链形式。

    四 资产组合模拟

    假设有五种资产,其日收益率(%)分别为

    0.0246 0.0189 0.0273 0.0141 0.0311

    标准差分别为

    0.9509 1.4259, 1.5227, 1.1062, 1.0877

    相关系数矩阵为

    1.0000 0.4403 0.4735 0.4334 0.6855

    0.4403 1.0000 0.7597 0.7809 0.4343

    0.4735 0.7597 1.0000 0.6978 0.4926

    0.4334 0.7809 0.6978 1.0000 0.4289

    0.6855 0.4343 0.4926 0.4289 1.0000

    假设初始价格都为100,模拟天数为504天,模拟线程为2,程序如下

    %run.m

    ExpReturn = [0.0246 0.0189 0.0273 0.0141 0.0311]/100; %期望收益

    Sigmas = [0.9509 1.4259, 1.5227, 1.1062, 1.0877]/100;%标准差

    Correlations = [1.0000 0.4403 0.4735 0.4334 0.6855

    0.4403 1.0000 0.7597 0.7809 0.4343

    0.4735 0.7597 1.0000 0.6978 0.4926

    0.4334 0.7809 0.6978 1.0000 0.4289

    0.6855 0.4343 0.4926 0.4289 1.0000

    ];%相关系数

    ExpCov = corr2cov(Sigmas, Correlations);%协方差

    StartPrice = 100;%初始价格

    NumObs = 504;

    NumSim = 2;

    RetIntervals = 1;

    NumAssets = 5;

    %开始模拟

    randn('state', 0);

    RetExact = portsim(ExpReturn, ExpCov, NumObs, RetIntervals,

    NumSim);

    Weights = ones(NumAssets, 1)/ NumAssets;

    PortRetExact = zeros(NumObs, NumSim);

    for i = 1:NumSim

    PortRetExact(:, i) = RetExact(:,:,i)*Weights;

    end

    PortExact = ret2tick(PortRetExact,

    repmat(StartPrice, 1, NumSim));

    plot(PortExact, '-r');

    蒙特卡洛模拟法简介

    蒙特卡洛(Monte

    Carlo)模拟是一种通过设定随机过程,反复生成时间序列,计算参数估计量和统计量,进而研究其分布特征的方法。具体的,当系统中各个单元的可靠性特征量已知,但系统的可靠性过于复杂,难以建立可靠性预计的精确数学模型或模型太复杂而不便应用时,可用随机模拟法近似计算出系统可靠性的预计值;随着模拟次数的增多,其预计精度也逐渐增高。由于涉及到时间序列的反复生成,蒙特卡洛模拟法是以高容量和高速度的计算机为前提条件的,因此只是在近些年才得到广泛推广。

    这个术语是二战时期美国物理学家Metropolis执行曼哈顿计划的过程中提出来的。

    蒙特卡洛模拟方法的原理是当问题或对象本身具有概率特征时,可以用计算机模拟的方法产生抽样结果,根据抽样计算统计量或者参数的值;随着模拟次数的增多,可以通过对各次统计量或参数的估计值求平均的方法得到稳定结论。

    二 蒙特卡洛模拟法求解步骤

    应用此方法求解工程技术问题可以分为两类:确定性问题和随机性问题。

    解题步骤如下:

    1.根据提出的问题构造一个简单、适用的概率模型或随机模型,使问题的解对应于该模型中随机变量的某些特征(如概率、均值和方差等),所构造的模型在主要特征参量方面要与实际问题或系统相一致

    2

    .根据模型中各个随机变量的分布,在计算机上产生随机数,实现一次模拟过程所需的足够数量的随机数。通常先产生均匀分布的随机数,然后生成服从某一分布的随机数,方可进行随机模拟试验。

    3.

    根据概率模型的特点和随机变量的分布特性,设计和选取合适的抽样方法,并对每个随机变量进行抽样(包括直接抽样、分层抽样、相关抽样、重要抽样等)。

    4.按照所建立的模型进行仿真试验、计算,求出问题的随机解。

    5. 统计分析模拟试验结果,给出问题的概率解以及解的精度估计。

    三 蒙特卡洛模拟法的应用领域

    蒙特卡洛模拟法的应用领域主要有:

    1.直接应用蒙特卡洛模拟:应用大规模的随机数列来模拟复杂系统,得到某些参数或重要指标。

    2.蒙特卡洛积分:利用随机数列计算积分,维数越高,积分效率越高。

    3.MCMC:这是直接应用蒙特卡洛模拟方法的推广,该方法中随机数的产生是采用的马尔科夫链形式。

    四 资产组合模拟

    假设有五种资产,其日收益率(%)分别为

    0.0246 0.0189 0.0273 0.0141 0.0311

    标准差分别为

    0.9509 1.4259, 1.5227, 1.1062, 1.0877

    相关系数矩阵为

    1.0000 0.4403 0.4735 0.4334 0.6855

    0.4403 1.0000 0.7597 0.7809 0.4343

    0.4735 0.7597 1.0000 0.6978 0.4926

    0.4334 0.7809 0.6978 1.0000 0.4289

    0.6855 0.4343 0.4926 0.4289 1.0000

    假设初始价格都为100,模拟天数为504天,模拟线程为2,程序如下

    %run.m

    ExpReturn = [0.0246 0.0189 0.0273 0.0141 0.0311]/100; %期望收益

    Sigmas = [0.9509 1.4259, 1.5227, 1.1062, 1.0877]/100;%标准差

    Correlations = [1.0000 0.4403 0.4735 0.4334 0.6855

    0.4403 1.0000 0.7597 0.7809 0.4343

    0.4735 0.7597 1.0000 0.6978 0.4926

    0.4334 0.7809 0.6978 1.0000 0.4289

    0.6855 0.4343 0.4926 0.4289 1.0000

    ];%相关系数

    ExpCov = corr2cov(Sigmas, Correlations);%协方差

    StartPrice = 100;%初始价格

    NumObs = 504;

    NumSim = 2;

    RetIntervals = 1;

    NumAssets = 5;

    %开始模拟

    randn('state', 0);

    RetExact = portsim(ExpReturn, ExpCov, NumObs, RetIntervals,

    NumSim);

    Weights = ones(NumAssets, 1)/ NumAssets;

    PortRetExact = zeros(NumObs, NumSim);

    for i = 1:NumSim

    PortRetExact(:, i) = RetExact(:,:,i)*Weights;

    end

    PortExact = ret2tick(PortRetExact,

    repmat(StartPrice, 1, NumSim));

    plot(PortExact, '-r');

    展开全文
  • 蒙特卡洛模拟

    2012-01-02 13:18:02
    蒙特卡洛(Monte Carlo)模拟是一种通过设定随机过程,反复生成时间序列,计算参数估计量和统计量,进而研究其分布特征的方法。具体的,当系统中各个单元的可靠性特征量已知,但系统的可靠性过于复杂,难以建立可靠性...
  • 通过考虑与测量误差和模型结构误差相关的不确定性,贝叶斯方法参数估计也可用于获得模型参数的置信区间。 然后,将蒙特卡洛模拟应用于模型参数的后验范围,以获取加氢裂化产物每个单独部分的模型输出的95%置信区间...
  • 蒙特卡洛法的基本思想是:为了求解问题,首先建立一个概率模型或随机过程,使它的参数或数字特征等于问题的解:然后通过对模型或过程的观察或抽样试验来计算这些参数或数字特征,最后给出所求解的近似值。...
  • 蒙特卡洛模拟法简介蒙特卡洛(Monte Carlo)模拟是一种通过设定随机过程,反复生成时间序列,计算参数估计量和统计量,进而研究其分布特征的方法。具体的,当系统中各个单元的可靠性特征量已知,但系统的可靠性过于...
  • 本文采用贝叶斯计算方法估计模糊数据的逆瑞利尺度参数。 基于不精确数据,无法以显式形式获得贝叶斯估计。... 此外,我们提供了通过蒙特卡洛模拟研究进行的数值比较,以均方误差值的形式获得比例参数估计值。
  • 上一节我们讲了蒙特卡洛方法,如果对蒙特卡洛方法不了解的朋友们可以看我之前写的这一篇文章。...序列估计问题谈起序列蒙特卡洛,就不得不说状态空间模型,可以说序列蒙特卡洛就是为了解决状态空间模型的参数估...
  • 蒙特卡洛法的基本思想是:为了求解问题,首先建立一个概率模型或随机过程,使它的参数或数字特征等于问题的解:然后通过对模型或过程的观察或抽样试验来计算这些参数或数字特征,最后给出所求解的近似值。...
  • 1、有向图与无向图模型 ...连续性:蒙特卡洛算法 3、生成模型和判别模型 CRF 属于判别模型 4、Log-linear model:一类模型 Logistic regression Conditional Random field 4.1 Multinomial Lo
  • Monte Carlo分析是一种器件参数变化分析,使用随机抽样估计来估算数学函数的计算的方法。它需要一个良好的随机数源。这种方法往往包含一些误差,但是随着随机抽取样本数量的增加,结果也会越来越精确。Monte Carlo ...
  • 马尔科夫链蒙特卡洛方法

    千次阅读 2013-10-14 19:15:47
    蒙特卡洛(Monte Carlo)模拟是一种通过设定随机过程,反复生成时间序列,计算参数估计量和统计量,进而研究其分布特征的方法。具体的,当系统中各个单元的可靠性特征量已知,但系统的可靠性过于复杂,难以建立可靠性...
  • 马尔科夫链蒙特卡洛(MCMC)

    万次阅读 多人点赞 2016-12-30 20:34:49
    在以贝叶斯方法为基础的机器学习技术中,通常需要计算后验概率,然后通过最大后验概率(MAP)等方法进行参数推断和决策。然而,在很多时候,后验分布的形式可能非常复杂,这个时候寻找其中的最大后验估计或者对后验...
  • 蒙特卡洛马尔可夫链学习小结

    千次阅读 2020-05-02 18:33:22
    MCMC背景 在统计学中,经常会遇到高维积分的计算,用传统的数值方法往往很难解决高维... 马尔可夫链蒙特卡洛MCMC,就是其中一种应用广泛的模拟方法,使得统计学领域中极大似然估计、非参数估计、贝叶斯统计学等...
  • 传染病模型的Matlab实现一、符号说明二、基本公式三、基本数据处理四、rβ参数估计及模型简化五、SIR模型MATLAB实现及分析六、蒙特卡洛方法的实现及分析七、模型改进八、文章说明 一、符号说明 以SIR模型为例,这里...
  • BES使用马尔可夫链蒙特卡洛(MCMC)估计贝叶斯框架中的参数。 BES包含一个通用模型,用于独立于交配系统来估计自交率和突变率。 BES还包含以下模型 纯雌雄同体 雌雄同体(雌雄同体+雄性) 妇科(雌雄同体+女性) ...
  • 蒙特卡洛模拟是一种通过设定随机过程,反复生成时间序列,计算参数估计量和统计量,进而研究其分布特征的方法。 蒙特卡洛方法的基本思想:所求解问题是某随机事件A出现的概率(或者是某随机变量B的期望值)。通过...
  • 蒙特卡洛(Monte Carlo)模拟是一种通过设定随机过程,反复生成时间序列,计算参数估计量和统计量,进而研究其分布特征的方法。具体的,当系统中各个单元的可靠性特征量已知,但系统的可靠性过于复杂,难以建立可靠性...
  • 本文是From Scratch: Bayesian Inference, Markov Chain Monte...马尔科夫链蒙特卡洛(MCMC, Markov Chain Monte Carlo)的定义是:通过在概率分布中进行采样,估计给定观测数据下模型的参数。(MCMC is a class of ...

空空如也

空空如也

1 2 3 4 5 ... 12
收藏数 235
精华内容 94
热门标签
关键字:

蒙特卡洛参数估计