• API
• FAQ
• Tools
• Archive
daily pastebin goal
63%
SHARE
TWEET

Modified python code

a guest Apr 4th, 2010 187 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. #!/usr/bin/env python
2.
3. from numpy.random import randn
4. from numpy import sqrt,row_stack,array,dot,transpose,sum,zeros,inner,zeros,mean,float32
5. def GMMGaussianMC(N=100):
6.     T=100;mu=0.1;sigma=0.2
7.     invT = 1./T
8.     coverage=zeros((2,N),dtype=int)
9.     d = array([[-1,0],[0,0]],dtype=float32)
10.     xmat = mu + sigma*randn(T,N)
11.     muhats =  invT * sum(xmat, axis=0)
12.     f1s = xmat - muhats
13.     sigma_hats = sqrt(invT * sum(f1s**2,axis=0))
14.     f2s = f1s**2 - sigma_hats**2
15.     f = row_stack((f1s,f2s))
16.     for n in xrange(0,N-1):
17.         #d is now inv(d) as it is easy to compute
18.         d[1,1] = (-2*sigma_hats[n])
19.         temp = f[:,n].reshape(2,T)
20.
21.         S=invT*inner(temp,temp)
22.         V=dot(dot(d,S),d)
23.
24.         SE_mu = sqrt(invT*V[0,0])
25.         SE_sigma = sqrt(invT*V[1,1])
26.         coverage[0,n]=(abs(muhats[n]-mu)<1.96*SE_mu)
27.         coverage[1,n]=(abs(sigma_hats[n]-sigma)<1.96*SE_sigma)
28.     y=mean(coverage,axis=1)
29.     return y
30.
31. if __name__=='__main__':
32.     #import cProfile
33.     #cProfile.run('GMMGaussianMC(100000)','prof_jason.dat')
34.     import time
35.     t1=time.time()
36.     y=GMMGaussianMC(100000)
37.     print time.time()-t1
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top