SHARE
TWEET

MLE Gaussian Example (fixed)

a guest Jan 28th, 2013 205 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import math
  2. import random
  3. from pylab import hist, show
  4.  
  5. #FUNCTIONS
  6.  
  7. def getNormalizerComponent(mu,sigma,dataRange):
  8. #make the resolution on this large.
  9.  n=0.
  10.  for count in range(dataRange[0]*100,dataRange[1]*100):
  11.   #density * volume
  12.   n+=(math.exp((-1.0/2.0)*pow(float(count)*0.01-mu,2)/pow(sigma,2)))*0.01
  13.  
  14.  return(n)
  15.  
  16.  
  17. #MAIN
  18.  
  19. print "Demonstration of MLE disentangling two mixed, truncated, unnormalized Gaussians"
  20.  
  21. #DECLARE THE TRUE DATA PARAMETERS (ANSWERS)
  22. trueSig1=0.9
  23. trueSig2=0.9
  24. trueMu1=4.
  25. trueMu2=8.
  26. trueRatio=0.2
  27. dataSize=100.
  28. range1=[-2,7]
  29. range2=[4,12]
  30.  
  31.  
  32. #DECLARE THE SEARCH SPACE
  33.  #brute grid
  34. sigma1List=[0.3,0.4,0.5,0.6,0.7,0.8,0.9]
  35. sigma2List=[0.3,0.4,0.5,0.6,0.7,0.8,0.9]
  36. mu1List=[1.,2.,3.,4.,5.,6.,7.,8.,9.,10.]
  37. mu2List=[1.,2.,3.,4.,5.,6.,7.,8.,9.,10.]
  38. ratioList=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9]
  39.  
  40. print "Generating Data..."
  41. #CREATE DATA ARRAY X
  42. x=[]
  43. #gauss1
  44. count=0
  45. while count < int(dataSize*trueRatio):
  46.  guess=random.gauss(trueMu1,trueSig1)
  47.  if guess > range1[0] and guess < range1[1]:
  48.   x.append(guess)
  49.   count+=1
  50.  
  51. #gauss2
  52. count=0
  53. while count < int(dataSize*(1.0-trueRatio)):
  54.  guess=random.gauss(trueMu2,trueSig2)
  55.  if guess > range2[0] and guess < range2[1]:
  56.   x.append(guess)
  57.   count+=1
  58.  
  59. #hist(x)
  60. #show()
  61. o=open("outfile.csv","w")
  62.  
  63. #LOOP THROUGH PARAMETERS IN A BRUTE FORCE GRID SEARCH
  64. maxProb=-99999999999
  65. for mu1 in mu1List:
  66.  print "Sanity Counter: "+str(mu1)
  67.  for mu2 in mu2List:
  68.   for sigma1 in sigma1List:
  69.    for sigma2 in sigma2List:
  70.     for ratio in ratioList:
  71.      workingProb=0.
  72.  
  73.      n1=getNormalizerComponent(mu1,sigma1,range1)
  74.      n2=getNormalizerComponent(mu2,sigma2,range2)
  75.        
  76.      for item in x:
  77.       prob1=0.; prob2=0.;
  78.       if range1[0] < item < range1[1]:
  79.        prob1=(ratio)*(1./n1)*math.exp((-1.0/2.0)*pow(item-mu1,2)/pow(sigma1,2))
  80.       if range2[0] < item < range2[1]:
  81.        prob2=(1.-ratio)*(1./n2)*math.exp((-1.0/2.0)*pow(item-mu2,2)/pow(sigma2,2))
  82.       workingProb+=math.log((prob1+prob2))
  83.  
  84.      o.write(str(mu1)+","+str(mu2)+","+str(sigma1)+","+str(sigma2)+","+str(ratio)+","+str(workingProb)+"\n")
  85.                  
  86.      if workingProb > maxProb:
  87.       maxProb = workingProb
  88.       maxSigma1 = sigma1
  89.       maxSigma2 = sigma2
  90.       maxMu1 = mu1
  91.       maxMu2 = mu2
  92.       maxRatio=ratio
  93.  
  94. o.close()
  95.  
  96. print "Values (true :: max)"
  97. print "Gaussian 1: ratio ("+str(trueRatio)+":"+str(maxRatio)+") // sigma ("+str(trueSig1)+":"+str(maxSigma1)+") // mu ("+str(trueMu1)+":"+str(maxMu1)+")."
  98. print "Gaussian 2: ratio ("+str(1.0-trueRatio)+":"+str(1.0-maxRatio)+") //sigma ("+str(trueSig2)+":"+str(maxSigma2)+") // mu ("+str(trueMu2)+":"+str(maxMu2)+")."
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. OK, I Understand
 
Top