欢迎使用CSDN-markdown编辑器

# -*- coding: utf-8 -*-"""Created on Mon May  8 10:20:25 2017@author: Administrator"""# -*- coding: utf-8 -*-"""Created on Mon Apr 24 15:50:10 2017@author: Dean"""#==============================================================================# There are 4 cols in array x, of which each col stands for k1, h1, k2, h2# All the variables are normal and the key parameters are:#    k1,5e-2,5e-3#    h1,1e-3,1e-4#    k2,4e-2,1e-3#    h2,15e-3,5e-4# Only care about h1 & h2#==============================================================================import numpy as npfrom scipy import stats as stsimport matplotlib.pyplot as pltplt.style.use('seaborn-colorblind')#==============================================================================##### MC Part#==============================================================================# Calculation module# Define the Response Surface functiondef RS(x):    xs = np.zeros(np.shape(x))    con1 = 0.7764490e+4*5e-2 - 0.3882245e+3    con2 = 0.7764490e+5*1e-3 - 0.7764490e+2    con3 = 0.7764490e+4*4e-2 - 0.3105796e+3    xs = 0.1552898e+4*x - 0.2329347e+2    y=0.2966269e+3+0.1892989e-2*con1-0.9514457e-1*con2+0.1998977*con3-0.5756422e+1*xs    y=y+0.2440628*xs**2+0.9268875e-2*con2*xs-0.8182638e-2*con3*xs    return y# Generate random samples for Monte Carlo simulations# Mean values & standard deviationmean = 15e-3std = 5e-4# Generate a sample populationns = 1000 # This number controls the smoothness of the sample populationsp = np.linspace(0.0001,0.9999,ns, endpoint = True)pop_mc = sts.norm.ppf(sp) * std + mean# Number of samples & experiments# Temperature tolerancetmp_ub = 310    # Upper boundary of temperaturent = 2000   # Calculation times in each repetitionrt = 1000   # Repetitions of simulationres_mc = np.zeros((nt,rt))intv = np.vectorize(np.int)sampleid_mc = np.round(np.random.uniform(0,ns-1,size = (nt,rt)))sampleid_mc = intv(sampleid_mc)sample_mc = pop_mc[sampleid_mc]res1 = RS(sample_mc)pos_mc = np.greater(res1, tmp_ub)   res_mc[pos_mc] = 1  # Select those exceeding tmp_ub and set as 1prob_mc = np.mean(res_mc, axis = 0) eff_mc = prob_mc    # Sampling efficiency, which equals to failure probability in MC sgm_mc = np.std(res_mc, axis = 0)   tmn_mc = np.mean(sgm_mc)/np.mean(prob_mc)   # Coefficient of variation#==============================================================================# ##### IS Part#==============================================================================opt_no = 20 # Optimization timestmn_trace = np.zeros(opt_no)mean_t = meanfor i in np.arange(opt_no):    dmean = np.random.uniform(-0.3,0) * std    mean_is = mean + dmean  # Change the mean value    pop_is = sts.norm.ppf(sp) * std + mean_is    sampleid_is = np.round(np.random.uniform(0,ns-1,size = (nt,rt)))    sampleid_is = intv(sampleid_is)    sample_is = pop_is[sampleid_is]    res2 = RS(sample_is)    res_is = np.zeros((nt,rt))    pos_is = np.greater(res2, tmp_ub)    res_is[pos_is] = 1    eff_is = np.mean(res_is, axis = 0)    alpha_u = sts.norm.pdf(sample_is[pos_is], loc = mean_t, scale = std)    alpha_l = sts.norm.pdf(sample_is[pos_is], loc = mean_is, scale = std)    alpha = alpha_u/alpha_l     # Get coefficient A_i    res_is[pos_is] = alpha    prob_is = np.mean(res_is, axis = 0)    sgm_is = np.std(res_is, axis = 0)    tmn_is = np.mean(sgm_is)    if(tmn_is<tmn_mc):  # Accept if better result acquired        mean = mean_is        tmn_mc = tmn_is        prob_best = prob_is    tmn_trace[i] = tmn_isplt.figure()plt.hist(prob_mc, 20, alpha = 0.8, label = 'MC')plt.hist(prob_best, 20, alpha = 0.8, label = 'IS')plt.legend()plt.grid('off')plt.title('Pf Estimation')plt.figure()#plt.hist(sgm_mc, 20, alpha = 0.8, label = 'MC')#plt.hist(sgm_is, 20, alpha = 0.8, label = 'IS')plt.hist(sgm_mc, 20, alpha = 0.8, label = 'MC', histtype = 'stepfilled')plt.hist(sgm_is, 20, alpha = 0.8, label = 'IS', histtype = 'stepfilled')plt.legend()plt.grid('off')plt.title('Std Comparison')plt.show()print(tmn_mc)print(mean)#==============================================================================