欢迎使用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 np
from scipy import stats as sts
import matplotlib.pyplot as plt
plt.style.use('seaborn-colorblind')
#==============================================================================
##### MC Part
#==============================================================================
# Calculation module
# Define the Response Surface function
def 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 deviation
mean = 15e-3
std = 5e-4
# Generate a sample population
ns = 1000 # This number controls the smoothness of the sample population
sp = np.linspace(0.0001,0.9999,ns, endpoint = True)
pop_mc = sts.norm.ppf(sp) * std + mean
# Number of samples & experiments
# Temperature tolerance
tmp_ub = 310    # Upper boundary of temperature

nt = 2000   # Calculation times in each repetition
rt = 1000   # Repetitions of simulation
res_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 1
prob_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 times
tmn_trace = np.zeros(opt_no)
mean_t = mean
for 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_is
plt.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)
#==============================================================================