Guest User

Modified python code

a guest
Apr 4th, 2010
223
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

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×