Demo entry 6356246

important sampling

   

Submitted by anonymous on Apr 17, 2017 at 17:55
Language: Python 3. Code size: 2.2 kB.

# -*- coding: utf-8 -*-
"""
Created on Wed Apr 12 14:29:10 2017
计算物理6 重要性取样
陈良玉    151120012
@author: think
"""

import numpy as np
import matplotlib.pyplot as plt

AA=0.0  #左边界
BB=20.0 #右边界
width=BB-AA  #区间长度
np.random.seed(456) #随机数种子

#所要求解的被积函数表达式
def ff(x):
    return x**4*np.exp(-x)
    
#权重函数1  
def weight1(x):
    return 1.0/width
    
#权重函数2
def weight2(x):
    return np.exp(-x)

#权重函数3
def weight3(x):
    return x**3*np.exp(-x)/6.0

#权重函数4
def weight4(x):
    return x**4*np.exp(-x)/24.0

maxSave=[1/width,1.0,0.3,0.3] #每个权重函数的最大值
#每个权重函数的引用
weightSave=[weight1,weight2,weight3,weight4]  

"""
下面对权重函数1和权重函数3情形进行多次计算
并给出结果的平均值,标准偏差,标准误差
"""
cal_num=25 #执行的总的计算次数
results1=[] #储存权重函数1情形的结果
results3=[] #储存权重函数3的情形的结果
results=[[],[]]
chooseIndex=[1,3] #即选择1和3两种情形
nt=2000  #每次计算时的样本容量

"""
下面计算产生results1,results3
"""
dd=0
while dd<np.size(chooseIndex):
    kk=chooseIndex[dd]
    mm=0
    while mm<cal_num:
        """
        mylist是用于储存所有的可以作为样本元素
        的delta的值,无需关心它们存入的先后顺序
        """
        mylist=[]
        i=0
        #产生符合指定分布的样本
        while i<nt:
            delta=AA+np.random.random()*width
            fvalue=weightSave[kk-1](delta)/maxSave[kk-1]
            randJudge=np.random.random()
            if randJudge<fvalue:
                mylist.append(delta)
            i+=1
         
        #按重要性取样计算积分值
        jj=0.0
        for ii in range(np.size(mylist)):
            x=mylist[ii]
            jj+=ff(x)/weightSave[kk-1](x)
        jj=jj/np.size(mylist)
        results[dd].append(jj) #储存结果
        mm+=1
    dd+=1

results1=results[0]
results3=results[1]

mean1=round(np.mean(results1),2)
std1=round(np.std(results1),2)
sterr1=round(std1/np.sqrt(cal_num),2)

mean3=round(np.mean(results3),2)
std3=round(np.std(results3),2)
sterr3=round(std3/np.sqrt(cal_num),2)

print('权重1对应的平均值',mean1,'权重3对应的平均值',mean3)
print('权重1对应的标准差',std1,'权重3对应的标准差',std3)
print('权重1对应的标准误',sterr1,'权重3对应的标准误',sterr3)

plt.plot(results1,'r',label='weight1')
plt.plot(results3,'b',label='weight3')
plt.ylim(20,28)
plt.xlabel('X')
plt.ylabel('Y')
plt.legend(loc='upper right')
plt.show()

This snippet took 0.00 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).