% matplotlib inline
from sympy import *
from scipy import random
import numpy as np
import matplotlib.pyplot as plt

a = 0
b = 1        
M = 100000   
xrand = np.zeros(M)
yrand = np.zeros(M)
zrand = np.zeros(M)

for i in range(len(xrand)):
    xrand[i] = random.uniform(a,b)
    yrand[i] = random.uniform(a,b)
    zrand[i] = random.uniform(a,b)

위에서처럼 코드를 짜서 xrand[i] = random value 형태로 i값에 따라 랜덤밸류가 달라지는 구조를 만들긴 했습니다만 그 이후가 생각보다 잘 되지 않습니다. 정확히 어떤 상황이냐면 다음과 같이.

import random

N = 900
def monte_carlo(the_list) :
    if len(the_list) > N :
        return the_list
    else :
        x = random.random()
        y = random.random()
        z = random.random()
        if p(y)/p(x) > z :
            the_list.append(y)    
            return monte_carlo(the_list)
        else :
            the_list.append(x)   
                # 재귀
            return monte_carlo(the_list)

new_list = []
print(monte_carlo(new_list)) 

위 코드는 예전에 제가 질문했을 때 어떤 고수분께서 올려주신 코드(나중에 알았지만 importance sampling 의 일종) 입니다.

이분은 재귀함수 형태로 뽑아주셨는데, 생각보다 오차가 크게 나와서 .. 방법을 좀 바꾸려고 합니다.

C로 코딩해서 성공한 연구실 동료가 조언해준 방법대로 i 값에 의존하는 x[i]를 만든 후 판정조건식 ( p(y[i])/p(x[i)] > z[i] ) 을 만족할 때는 그대로 y를 append 하고,

else의 경우에는 재귀를 취하지 않고, 그냥 i = i + 1 로 넘어가서 ( i + 1)에서 새로운 판정 p(y[i+1]) /p(x[i+1]) > z[i+1] 을 해나아가는.. 그런 구조를 만들어서, 판정을 만족시킨 N개의 random value들을 얻고 싶은데 잘 되지 않네요.. (참고로 이 경우에도 M개의 random value 값들이 N개의 random 값들을 추출하기 전에 소진될 가능성이 있습니다. 무한개의 random value 리스트를 만드는게 가능한지 몰라서 일단은 이 방법으로 해보려고 합니다.)

어떻게 해야 할까요?... 조언 부탁드립니다.

추가) 일단 제가 며칠전에 완성했던 허졉한 코드는 다음과 같습니다.. 부끄럽지만 조언 부탁드려여 ㅠ.ㅠ 코드의 중간 부분부터가 importance 샘플링입니다.

from sympy import *
x = Symbol('x')
y = Symbol('y')
z = Symbol('z')
a = Symbol('a')

def f(x):
    return 3 * x **2
def p(x):
    return 2 * x 
def g(x) :
    return f(x)/p(x)

Integral(f(x),x).doit()
pprint(Integral(f(x),x).doit()) 
pprint(Integral(f(x), (x,0,1)).doit()) # accurate integration value
pprint(Integral(p(x), (x,0,1)).doit()) # Normalization check

# 구분구적법으로 함수 정적분 풀이 시도

% matplotlib inline
from sympy import *
from scipy import random
import numpy as np
import matplotlib.pyplot as plt

a = 0
b = 1        # limit of integration 0부터 파이까지 함수 적분
M = 100000   # 쪼개는 수
xrand = np.zeros(M)

for i in range(len(xrand)):
    xrand[i] = random.uniform(a,b)

integral = 0.0
for i in range(M):
    integral += f(xrand[i])

answer = (b-a)/float(M)*integral
error = abs(1-answer)
print("The integral from 0 to 1 of f(x): ", answer)
print("error: ", error)


import random

# importance sampling 시도해본다.

'''
 엽토님께.

 기본적인 아이디어는 integral(f(x))dx(0부터 1까지) 라는 정적분값을 난수를 이용해 근사과정으로 구하는 것입니다.
 f(x) = (f(x)/p(x)) * p(x) = g(x) * p(x) 라 하고 이때 p(x)를 weight function 또는 확률분포함수(pdf)라고 합니다.

 이 때 특정한 판정조건(아래코드의 조건식)을 만족하는 x값들을 N개 추출하고 이를 통해 g(x)의 평균값을 구하면 
 " Mean value of g(x) = f(x)의 원래 정적분값 " 이 성립함이 알려져 있습니다. N값이 커질수록 두 값의 차이는 줄어듭니다.

 (참고로 이 문제에서는 f(x) = 3x^2, p(x) = 2x 라 가정하여 g(x) = 3x^2/2x로 되어 있다고 보시면 됩니다.)
 (따라서 최종적으로 가까워져야 할 값은 정적분값인 1입니다.)
''' 

# 비어 있는 리스트를 주면 특정 식을 만족할 경우 해당 랜덤값을 최대 N개까지 목록에 넣어서 돌려준다.
N = 900
def monte_carlo(the_list) :
    if len(the_list) > N :
        return the_list
    else :
        x = random.random()
        y = random.random()
        z = random.random()
        if p(y)/p(x) > z :        # p(X)는 이 경우 1차 다항함수이므로 y/x 해도 상관없음.
            the_list.append(y)    # 이 함수 자체를 다시 실행
            return monte_carlo(the_list)
        else :
            the_list.append(x)   # x를 리스트에 꼽고
                # 재귀
            return monte_carlo(the_list)

new_list = []
#print(monte_carlo(new_list)) 

I = 0
for i in range(N):
    I += g(monte_carlo(new_list)[i])/ N
print(I)